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Abstract. Monotone piecewise cubic interpolants are simple and effective. They 
are generally third-order accurate, except near strict local extrema where accuracy de- 
generates to second-order due to the monotonicity constraint. Algorithms for piecewise 
cubic interpolants, which preserve monotonicity as well as uniform third and fourth- 
order accuracy, are presented. The gain of accuracy is obtained by relaxing the mono- 
tonicity constraint in a geometric framework in which the median function plays a 
crucial role. 

1. Introduction. When the data are monotone, it is often desirable that the 
piecewise cubic Hermite interpolants are also monotone. Such cubics can be constructed 
by imposing a monotonicity constraint on the derivatives, see [2], [3], [5]— [10] , [17], [19]. 
These constraints are obtained from sufficient conditions for a cubic to be monotone. 
A necessary and sufficient condition for monotonicity of a Hermite cubic was indepen- 
dently found by Fritsch and Carlson [10], and Ferguson and Miller [8]. The resulting 
constraints are complicated and nonlocal [7], [10]. A simpler sufficient condition, which 
generalizes an observation by De Boor and Swartz [5], is more commonly used [10], [17], 
[19]. The corresponding constraint is local. It imposes different limits on the derivatives 
depending on whether the data are increasing, decreasing, or have an extremum. In 
[17], Hyman simplified De Boor and Swartz’s method by enforcing a single limit for all 
cases. Hyman’s constraint is also less restrictive and slightly more accurate. However, 
in some special cases as shown in subsection 2.3 below, this constraint may create os- 
cillations. A different idea for simplification was presented by Butland [2], and Fritsch 


and Butland [9]. They constructed limiter functions whose arguments are the slopes of 
the piecewise linear interpolation. The resulting derivatives satisfy the monotonicity 
constraint a priori , but they are at best second-order accurate. The constraints of De 
Boor and Swartz and Hyman, on the other hand, can be applied to any high-order 
approximation of the derivatives. 

A major problem of the above methods, including the fourth-order algorithm in [7], 
is that if the data are no longer monotone, the constraints may cause a loss of accuracy. 
In fact, near strict local extrema, the interpolants are only second-order accurate, i.e., 
the error is 0(/i 2 ), and thus, no better than linear interpolation. Recently, Hyman’s 
constraint was modified to obtain uniform third-order accuracy [6]. The proofs of 
stability and monotonicity, however, were not shown. 

In a different context, namely numerical simulations of conservation laws, there 
has been a great deal of effort to preserve monotonicity by constraining the deriva- 
tives, see [12], [18], [20]— [22], [24], [27], and the references given there. In an earlier 
paper [26], Van Leer introduced a monotonicity constraint which corresponds to De 
Boor and Swartz’s constraint for quadratic interpolation. His harmonic mean limiter 
is also identical to that of Butland [2]. In this context, the study of limiter func- 
tions was simplified by reducing them to functions of one variable instead of two. It 
is well-known that all limiters suffer a loss of accuracy near strict local extrema. Re- 
cently, Roe [22] suggested that “Although only very simple ideas are involved, more 
research is still needed to clarify the properties of limiter functions.” In [15], Harten 
and Osher solved the loss of accuracy problem by combining monotonicity and accu- 
racy into nonoscillatory, piecewise parabolic interpolation. The concept behind their 
UNO method completely deviates from that of the monotonicity constraint: it defines 
high-order accurate derivatives which do not cause oscillations, but it does not provide 
a constraint for other approximations of the derivatives. The ENO schemes [14], which 
are higher-order extensions of the UNO method, are unstable in the sense that a small 
change in the data can produce a large change in the interpolant. Although the UNO 
and ENO schemes are highly accurate, they are somewhat diffusive, e.g., they smear 
contact discontinuities [13], [28], [30]. Modifications of these schemes to improve the 
resolution at contact discontinuities can be found in [13], [28]. In [16], the present 
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author took a different approach to obtain uniform high accuracy: the monotonicity 
constraint was extended in a geometric framework using the UNO method. The re- 
sulting SONIC schemes axe highly accurate [16], [29], and the UNO scheme itself is a 
member of the SONIC class. 

We follow the author’s geometric approach below. First, the concept of monotonic- 
ity of the data in an interval is introduced so that monotone parts can be distinguished 
from strict local extrema. Next, we present a geometric framework in which the median 
function plays a crucial role for all methods. In this framework, De Boor and Swartz’s 
constraint is essentially as simple as Hyman’s. Since the latter may cause oscillations, 
we choose to extend the former. Our systematic approach results in several uniform 
third-order as well as fourth-order constraints which are stable and can be applied to 
any high-order approximation of the derivatives. The geometric framework and the 
introduction of new concepts also made possible the proofs of monotonicity, accuracy, 
and stability. We note that the UNO scheme is a member of the class of third-order 
methods presented here. It produces interpolants which are the “flattest” in this class. 
On the contrary, the fourth-order ENO scheme, which is unstable, does not belong to 
the fourth-order class whose members are stable. Since several of our methods employ 
limiter functions, we present a rather complete analysis of the essential properties of 
limiters. Finally, for smoother transitions near the boundaries, we introduce nonlinear 
boundary conditions via the median function. Note that our methods can also be com- 
bined with cubic splines in a manner similar to Hyman’s method [17]. However, only 
local methods are presented below. 

This paper is essentially self-contained. In section 2, monotonicity of the data in an 
interval is defined, and well-known results on Hermite cubics are reviewed. The De Boor 
and Swartz monotonicity constraint and algorithm are presented using our framework 
in subsection 2.1. Subsection 2.2 deals with limiter functions by first reducing them 
to functions of one variable. It is shown that the construction of a limiter g can be 
reduced to that of a function on [0,1] with #(0) = 0, and < 7 ( 1 ) = !• Furthermore, 
the order of accuracy of a limiter is determined by the slope <7 f (l) in the case of a 
uniform mesh. In subsection 2.3, Hyman’s method is described, and examples which 
show the difficulties of obtaining both accuracy and monotonicity axe pxesented. These 
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difficulties are then resolved in section 3. Subsection 3.1 describes the UNO scheme 
using the median concept. In subsequent subsections, three different monotonicity- 
preserving constraints, which are stable and provide third-order accurate interpolants, 
are presented. It is shown that these constraints reduce to that of de Boor and Swartz 
for monotone data; however, they can be less restrictive when the data are no longer 
monotone. In subsection 3.5, besides additional analysis and comparison of third- 
order methods, we present examples which show that all third-order constraints may 
not reproduce a cubic, and that near inflection points with zero slopes, they are no 
higher than third-order accurate. In section 4, fourth-order methods are presented. 
Nonoscillatory cubics, which are introduced in subsection 4.1, are employed in the 
definition and the development of fourth-order q-monotonicity-preserving constraints 
in subsection 4.2. Several nonlinear boundary conditions are presented in section 5 
by using the median function. Numerical results axe shown in section 6. Finally, the 
discussion on some applications of the new methods and the conclusions are presented 
in section 7. For readers who are only interested in high-order monotonicity-preserving 
constraints, we suggest that they skip subsections 2.2 and 2.3. 

2. Cubic Hermite interpolation. Let the mesh {x,}"_ m , m < n, be a partition 
of the interval [x m ,x n ] such that x m < x m +i < ... < x n _i < x„. Let {fi} be 
the corresponding data which are samples of a piecewise smooth function /, and let 
Fi = (x;, fi ). The local mesh spacing is Ax l+1 / 2 = x l+ i — x i, and the maximum value 
of all Axj +1 / 2 is denoted by h. The slope of the piecewise linear interpolant is equal to 
the first divided difference, 

s i+l/2 — /[*,-, x j+1 ] — A/, + i/ 2 /Ax, + 1 / 2 . (2-1) 

The data are increasing at Xj if /;_i < /, < /<+ j. They are increasing in [x;,x,+i] if 
they are increasing at X{ and Xj+j, i.e., /j_i < /, < fi+% < fi+ 2 - Similar definitions 
hold for decreasing data with appropriate sign changes. The data are monotone at 
Xi (or in [xj,x 1+ i]) if they are increasing or decreasing at x, (or in [xj,x, + i]). Note 
that our definition of monotonicity of the data in [xj,x;+i] is rather straightforward; 
however, it plays an essential role in the high-order extensions of the monotonicity 
constraints and their proofs. The interpolant Pf is monotone in [x;,x, +1 ] if (P/)(x) 
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is monotone for x between Xi and i,-+ j. The interpolant is of class C k if (Pf)(x) is 
continuous and has continuous derivatives for all orders less than or equal to k. 

Given the data {/,} and the slopes {/J, which respectively approximate the exact 
values {/(x;)} and {/'(xj)}, the cubic Hermite interpolant is defined for x m < x < x n 
in terms of fi+i, fi, fi+i by (see [4] or any standard text on numerical methods), 

q(x) = c 3 (x - x,) 3 + c 2 (x - x,) 2 + ci(x — x<) + cq (2.2) 

where, for x, < x < x l+1 , 

co — fi > Cj fi , 

c 2 = (3^, -1-1/2 — 2/j — /i+l)/(AXi_j.|/ 2 ), 

C3 = (fi + fi+1 ~ 2Si + i/ 2 )/(Axi + i/ 2 ) 2 . 

One can easily verify that q(xi) = fi, q'(xi) = fi, q(x i+1 ) = / I+1 , g'(x,- + i) = f i+ j. 
Consequently, the interpolant (2.2) has a continuous first derivative (<?(x) € C 1 ). It is 
fourth-order accurate if the derivatives {fi} are exact or third-order, and third-order if 
the derivatives are second-order, etc. These derivatives can be approximated by using 
polynomial interpolation. If quadratic interpolation is used, one obtains the following 
formula which is second-order accurate: it approximates the exact derivative to 0(h 2 ), 

j. _ A J »~l/2 3 »+l/2 + Ax t+ i / 2 3,_ 1 / 2 

Ax,_i/ 2 + Ax i+ i/ 2 

(Note that the order of accuracy as defined in the literature dealing with conservation 
laws is one less than that defined here, which follows the definition in the literature 
dealing with monotone interpolation). The above expression is called the parabolic 
formula [17] since it gives the slope at x, of the parabola P,(x) = (x,p,(x)) through the 
points Fi- 1 , Fi, and F,+i. Such an approximation has the advantage of being local in 
the sense that only nearby data are used. The resulting interpolant (2.2) is also local. 
If the data {fi} are monotone, the interpolant is not necessarily monotone. One needs 
additional constraints on {/,}. 

2.1. Monotonicity constraints. Let I[z \, . . . , x*] be the smallest closed interval 
containing Z \ , . . . , 2*, i.e., 

I [ z i , ■ • • , z k\ — [min(z 1 ,...,z*),max(x 1 ,...,xj k )]. 
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Let the median of three numbers be the one which lies between the other two. Notice 
that the median lies in the interval defined by any two of the three arguments. Let 


minmod (x, y) = median (x, y, 0), 


then a convenient formula for the minmod function is 


minmod (x,y) = ~ [sgn(x) + sgn(y)] min(|x|, |y|) 


(2.4) 


where sgn(x) = 1 if x is positive, and sgn(x) = -1 if x is negative. Observe that if 
x = 0, it does not matter whether sgn(x) is defined as -1, 0, or 1 in (2.4) since the 
resulting minmod is 0. The minmod function is more commonly defined as, see e.g., 

[12], [15], [24], [26], 

minmod (a:, y) = { “™(M.lvl) if *V >0, 

1 0 otherwise. 


Conversely, the median function can be expressed in terms of minmod: 

median (x, y, z) = x + minmod (y — x, z — x) 

= y + minmod (x — y, z — y). 


(2.5) 


The following sufficient condition for monotonicity is a generalization of an observation 
by De Boor and Swartz [5]. It was developed and used by many authors [6], [10], [17], 

[19]- 


Lemma 1. If 

fi, fi+i 6 /[0,3s i+1 / 2 ]> (2-6) 

then the resulting interpolant (2.2) is monotone in [xj,x i+1 ]. 

The constructive proof below is simpler than those of the necessary and sufficient 
condition [6], [10]. The case of fi = fi+\ is trivial since the interpolant is a constant 
function. If fi fi+i, consider the linear change of coordinates x = (x — x ( )/Ax, +1 / 2) 
y = (y — fi)/ Afi+i/i, where (x,y) is the original coordinate system. Clearly, the cubic 
y = q(x) defined by (2.2) is transformed to a cubic y = q(x). By using the chain rule, 
the monotonicity of q(x) for x, < x < x 1+J is equivalent to the monotonicity of q(x) 
for 0 < x < 1. Without loss of generality, therefore, one may assume x, = 0, x I+1 = 1, 
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and fi = 0, /,+ 1 =1. As a result, each pair (/,-,/,+i) determines an interpolant. The 
four pairs (fi,fi+i) = (0,0), (3,0), (0,3), and (3,3) form the comers of a square called 
the De Boor and Swartz square shown in Fig. la. Denote the four corresponding 
interpolants by /q 0 ,o)> ^( 3 ,o)> fyo, 3 )> and ^(3,3)- Using definition (2.2), it follows that 


h( 0 ,o)(x) = -2x 3 + 3x l 2 , h( 3 ,o)(x) = x 3 - 3x 2 + 3x, 
V, 3)0*0 = x 3 , h (3)3 )(x) = 4x 3 - 6x 2 + 3x. 



(a) De Boot and Swartz’s square (b) The interpolants of its four corners 
Figs. 1. De Boor and Swartz’s square and the interpolants of its four corners. 


One can verify that the above cubics are increasing for 0 < x < 1, see Fig. lb. Any 
linear combination of them with positive coefficients is also increasing for 0 < x < 1. 
Denote (/,, /,+i) by (a,/3), then condition (2.6) is equivalent to 0 < a/3 < 1, and 
0 < /?/3 < 1. Consequently, the following cubic Hermite interpolants are increasing for 
0 < x < 1: 

, a a. 

«(a,0) ~ ■g *(3,0) + (1 - )h(0,0), 

l a a 

"(or, 3) = 2 "(3,3) + (1 ~ -g )«(0,3)> 

h(a,p) = f fya,3) + (1 - ^)^(a.O)- 

Note that the slopes (/,-,/,+i) of the above cubics are respectively (a, 0), (a, 3), and 
(a, /?). This completes the proof. 
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Ferguson and Miller [8], and Fritsch and Carlson [10] independently found a neces- 
sary and sufficient condition which is an extension of (2.6), see also [6]. Because of its 
simplicity, however, condition (2.6) is more commonly used. Note that in the context of 
conservation laws, the factor 3 in (2.6) is replaced by 2. The reason for this difference 
is that the solutions, as functions of time steps, correspond to quadratic rather than 
cubic interpolations, see [16]. 

Expression (2.6) with index i replaced by i — 1 together with (2.6) imply 

fi G -f [0, 3sj_x/2] and /, G /[0, 3sj+i/2]- (2-6 ) 

Consequently, /,• belongs to the intersection of /[0, 3 s,_i/ 2] and / [0, 3s t+1 /2], which is 
/ [0, minmod (3s,_ i/2,3sj +1/2 )]. Thus, if 

/, G I [0, 3s,-] (2.7) 

where 

Si = minmod (*,-- 1 / 2 . 01 + 1 / 2 ) (2.8) 

for all i , m < i < n, then the interpolant (2.2) is monotone in [x,,x,+i], m < i < n - 1. 
To bring /, into the interval I [0, 3s,], it suffices to replace /, by the median of /,, 

0, and 3sj. Using the definition of the minmod function, one obtains the monotonicity 
algorithm 

/i <— minmod (/i,3si). (2-9) 

Notice that (2.9) can also be expressed as 

fi <— minmod[minmod(/,,3sj_ 1 /2),3s,- + i/2], (2-9') 

1. e., the constraints are enforced on the left, and then on the right. Both (2.9) and 
(2.9') are equivalent to the following commonly used formula: 

{ 0 if •S,-l/2 5 «+l/2 ^ 0) 

min[max(0,/i),3min(s,_i/ 2 ,5,+i/2)] if -S.- 1/2 and s i+1 / 2 > 0, 
max[min(0,/,),3max(s,_i/2,s i+1 / 2 )] if 1/2 and s i+ i/ 2 < 0. 

Expressions (2.9), (2.7), and the interval in (2.7) are called the MP (monotonicity- 
preserving) algorithm, the MP constraint, and the MP interval, respectively. The 
combination of (2.3) and (2.9) is called the MP-parabolic algorithm. 


8 



Next, we show that if / € C 2 , {fi} are second or higher-order accurate, and 
the original {,/i} are at least first-order, then at the part where the exact derivative 
f'(x) f 0, the MP algorithm has no effect provided the mesh is fine enough. Indeed, let 
xi be in a small neighborhood of x such that f'(xi) ^ 0. Since both s,-i/ 2 and 
approximate f'(xi) ^ 0 to 0(h), they are both strictly positive, or strictly negative for 
small enough h\ consequently, the MP interval contains /[0, 3/ ; (xi)] up to an error of 
the order 0(h). Because the original f, approximates f'(xi) to at least 0(h ) , it belongs 
to the MP interval, and thus is not altered by the MP algorithm. 

To analyze the MP constraint near strict local extrema, assume that / £ C , {/,} 
are third or higher-order accurate, and the original {fi} are at least second-order. Near 
a strict local extremum ati = i*, suppose f"(x*) / 0. The following arguments hold 
in a small neighborhood of x* where higher-order terms can be neglected. Again, by 
using a linear change of coordinates as shown in the proof of Lemma 1, we may assume 
x* = 0, and f(x) = x 2 + higher-order terms. We claim that if 0 < x,_i < x,- < x,+i 
(similarly for x j <C x,' < x j ^ 0), then the MP algorithm does not alter fi provided 
h is small enough. Indeed, by calculating £j-i/ 2 > ^t+i /2 from the above expression for 
/(x), it follows that the MP interval contains [0, 3x*]. Since f'(xi) = 2 x, + 0(x?), and 
\f i — / f (xj)| = 0(x 2 ), one concludes that fi lies in the MP interval if h is small enough, 
and the claim follows. Observe that there are at most two data points adjacent to the 
local extremum which satisfy neither x* < Xj-i < Xi < x,-+i nor Xi_i < Xi < x,-+i < x . 
At these points , x,_j < x* < Xj+i, and the MP algorithm may change the original 
fi, causing a loss of accuracy. 

Note that in the MP-parabolic algorithm, both (2.3) and (2.9) are defined in terms 
of s;_i /2 and s i+1 / 2 - Therefore, one may attempt to derive formulas for fi which satisfy 
the MP constraint (2.7) a priori. Such an approach is presented below. 

2.2. Limiter functions. For convenience of notation, in the rest of this section, 
denote Sj_i/ 2 and Sj +1 y 2 by s and t , respectively. Let 

fi = G(s,t) 

where G is a nonlinear averaging function called a limiter function. Fritsch and Butland 
[9] analyzed G as a function of two variables. We follow the analysis commonly used 
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in conservation laws by first reducing G to a function of one variable, see e.g., [18], 
[20]— [22] , [24], [26], [27]. The analysis below is more complete and leads to many new 
limiters. In order for G to give a result independent of units of measurement, it is 
necessary that G is homogeneous of degree one, that is, for any real number k, 

G(ks,kt ) = kG(s,t). 

Substitute k = 1/t into the above equation, one obtains 

G(s, t ) = tG(s/t, 1) = iG(r, 1) = tg(r) (2.10) 

where r = s/t and g(r ) = <7(r, 1). Expression (2.10) shows that G is completely 
determined by g ; therefore, g is also called a limiter. Since G is an averaging function, 
the following properties for G are desirable (see also [9]): 

a. G is symmetric, i.e., (7(M) = G(t,s). FYom (2.10), this is equivalent to 

tg(r) = sg(l/r), or </(l/r) = (l/r)g(r). (2.11) 

b. (j(.s, t) lies between s and t. Equivalently, 

0 (r)€/[l,r], (2.12) 

i.e., the graph of g lies in the shaded region of Fig. 2a, which is called the 
second-order region (more precisely, sufficient for second-order accuracy re- 
gion) for the following reason: since both s and t approximate /'(a;, ) to 0(h ) , 
the above condition implies the same is true for /*; consequently, the in- 
terpolant (2.2) is at least second-order accurate. Note that (2.12) implies 
<7(1) = 1- 

The next two conditions correspond to stability and monotonicity: 

c. G or g is continuous, 

d. C7(s,t) satisfies the MP constraint (2.7), i.e., 

g(r) € I[0,3minmod(l,r)]. (2.13) 

Observe that the above condition implies that if r < 0, then g(r) = 0. For 
r > 0, it means that the graph of g lies in the shaded region of Fig. 2b. This 
region, together with the negative r-axis, is called the MP region. 
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(c) MP2 region (d) Average limiter 

Figs. 2. Limiters. 

In this paper, unless otherwise stated, all limiters satisfy the first three of the 
above four conditions, and limiters in this subsection also satisfy the fourth. Note that 
in the context of conservation laws, limiters need not satisfy all of the above conditions, 
see e.g., [18], [23], [24]. 

For r > 0, condition (2.12) together with (2.13) axe equivalent to the fact that the 
graph of g lies in the intersection of the shaded regions of Figs. 2a and 2b, which is 
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shown in Fig. 2c. This region, together with the negative r-axis, is called the MP2 
region. (The MP2 region is similar to the TVD2 region defined by Sweby [24]; the 
difference is that the constant 2 is replaced by 3.) Furthermore, for r > 1, g(r) can 
be defined in terms of g(r ) where 0 < r < 1 thanks to the symmetric property (2.11). 
Thus, the task of defining the limiter function G reduces to that of defining a continuous 
function g(r) for 0 < r < 1 such that <7(0) = 0, <7(1) = 1, and the graph of g lies in the 
region bounded by the three lines <f>(r ) = r, <f>(r) = 1, and <f>(r) = 3 r. Some popular 
limiter functions are listed below. 

The minmod limiter [12], which corresponds to the lower boundary of the MP2 
region, is defined by 

g(r) = minmod (l,r), or G(s,t) = minmod (s, t). (2-14) 


To obtain the harmonic mean limiter, which was discovered independently by Van 
Leer [vL] and Butland [2], one represents g by a rational function of the form ar/(br + 1) 
where a and b are constants. If we assume g'(0 + ) = 2, then a = 2; since <7(1) = 1, 
6=1. Thus, 


ff( r ) = 


r + |r| 

1 + M’ 


(2.15a) 


or 


G(s,t) 


' 0 

< 2 st 

< i 5 + < 


if st < 0, 
if st > 0. 


(2.156) 


Fritsch and Butland [9] observed that the harmonic mean limiter produces curves 
that are “too flat” because <?'(0 + ) = 2. They proposed the following limiter function 


with 0'( o+ ) = 3, 


9(r) 


' 0 

3r 

< l + 2r 
3r 

k 2 + r 


if r < 0, 
if 0 < r < 1, 

if r > 1, 


(2.16) 


or 


G(s,t) 


’ 0 

3s t 

* 2s -)■ t 
Zst 

. s -J- 2 1 


if st < 0, 

if st > 0 and |s| < |t|, 
if st > 0 and |s| > |t|. 
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The limiter which corresponds to the upper boundary of the MP2 region is similar 
to Roe’s “Superbee” limiter [20], [24]: 


g(r) = min{max(l,r), 3 minmod (l,r)}. 


(2.17a) 


The corresponding expression for G(s,t ) is obtained by first defining the maxmod 
function, 

maxmod ($,t) = - [sgn(s) + sgn(<)] max(|s|, \t\), 

G(s,t) = minmod [maxmod (s,<),3minmod(s,<)] 

1 . (2.175) 

= - [sgn (s) + sgn(f)] min [max(|s|, |t|),3mm(|s|, |<|)]. 

The above definition of the maxmod function is unstable when s = 0 or t = 0; however, 
since the limiter returns 0 in these cases, it is stable. Note that our definition of 
the maxmod function simplifies the expression of the “Superbee” limiter; furthermore, 
(2.17b) extends to higher-order cases in the next sections. 

The minmod, the Fritsch-Butland, and the “Superbee” limiters are only first-order 
accurate because the slopes of these limiters are discontinuous at r = 1 (more on this 
later). Assume that g is piecewise C 1 , then g' is continuous at r = 1 if and only if 
g\ 1") = 1/2, or g'{ 1+) = 1/2, or 

ff'(l) = 1/2. (2-18) 

Indeed, by differentiating (2.11), 

s ' (r) + i s ' (i) = s (;) . 

and the above claim follows, since g(l) = 1. Furthermore, if g is piecewise C 2 , then 
again by differentiating the above equation, 

© ■ 

Consequently, ^"(l - ) = g"{ l + )> i-e., g" is continuous at r = 1 a priori. We present 
below three limiter functions with (?'(1) = 1/2 and 5f'(0 + ) = 3. 

The average limiter is defined by enforcing the constraint on the average, (see [26], 
[27] for the average limiter with <?'(0 + ) = 2), 


g(r) = minmod | — ■t — ,3 minmod (l,r)|^ , 


(2.19) 
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or 


G(s,t) = minmod 



3 minmo 



The above limiter is the combination of Fromm’s scheme [11] and the MP algorithm. 
Its graph is shown in Fig. 2d. 

If we use a rational function of the form (o^r 2 + ajr)/(r 2 + br + 1) to calculate 
g , then, in order that </(()) = 3, <zi = 3. Since g tends to 3 as r tends to oo, 02 = 3. 
Finally, </(l) = 1 implies 6 = 4, 


or 


g(r) = 


0 

3r 2 + 3r 
r 2 + 4r + 1 


if r < 0, 

if r > 0, 


G(s,t) 


0 

3 st(s + t ) 
s 2 + 4 st + t 2 


if st < 0, 
if st > 0. 


( 2 . 20 ) 


If a cubic of the form (2.2) is used, the result is 


9{r) 


0 

3 3 7 2 „ 

r ~2 r +3r 

6r 2 - 7r + 3 
2r 2 


if r < 0, 
if 0 < r < 1, 

if r > 1. 


( 2 . 21 ) 


Since the corresponding expression for G($,t) can be obtained in a similar manner as 
the above limiters, it is omitted. 

Next, we show that the harmonic mean and the above three limiters are second- 
order accurate. 

Lemma 2. Suppose / € C 3 , and the data {/, } are third or higher-order accurate; the 
mesh is uniform; the limiter function g is continuous, piecewise C 2 , and satisfies (2.11) 
and (2.12); then at the part where f'(x) / 0 , /, defined by g approximates the true 
derivative /'(x, ) to second-order accuracy for small enough h if and only if g'( 1) = 1/2. 

First, suppose y'(l) = 1/2. We will show that f, is second-order. Since the mesh 
is uniform, the limiter g a (r) = (r + l)/2 yields identical /; as the parabolic formula 
(2.3); consequently, it is second-order. Because both s and t approximate /'(x* ) ^ 0 
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to 0(h ) , if h is small enough, t is bounded away from 0, and r - 1 = 0(h). Since g a 
and g have the same slopes at r = 1, 

g(r) - g a (r) = 0((r - l) 2 ) = 0(h 2 )-, 

therefore, /, defined by g is second-order. 

Next, suppose = 0 / 1/2. We will provide an example where the resulting 

fi is only first-order accurate (the case </(l + ) ^ 1/2 is similar). Let /, be defined by 
the quadratic f(x) = x 2 -f x, and let x,- = 0, i = —h, and x*+i = h. Then s = 1 — h, 
t = l + h, r = l — 2h + 0(h 2 ). From (2.10), 

fi = (1 + h)[ 1 + 0(r - 1) + 0(h 2 )] = 1 + (1 - 2 0)h + 0(/i 2 ). 

Since /'(x,) = 1 and 0 / 1/2, the above /, is first-order. This completes the proof. 

For an irregular mesh, all limiters are first-order accurate, while the MP-parabolic 
algorithm is still second-order away from strict local extrema. Near strict local extrema, 
however, the monotonicity constraint causes a loss of accuracy as shown below. 

2.3. Accuracy and monotonicity. If the data are not monotone in [x,-,x,+i], 
the MP algorithm (2.9) still produces a monotone interpolant, which causes a loss of 
accuracy as in the following well-known example, see e.g., [6]. Let the data points be 
on a parabola such that /, = / J+ i as shown in Fig. 3; then the MP constraint (2.7) 
implies that fi = /,+ 1 = 0. Therefore, /, and /,+i are only first-order accurate, and 
the cubic interpolant reduces to linear interpolation between fi and /<+ 1 - Observe that 
the interpolant is still monotone in [x,, Xj+i] even though the data are not. One may 
attempt to turn off the monotonicity constraint near strict local extrema to avoid the 
above “clipping” phenomenon; however, this makes the method unstable in the sense 
that a small change in the data may cause a large change in the interpolant. 

In [17], Hyman extended the MP constraint to 

fi <E /[— 3min(]s|, |#|),3min(|s|, |t|)]. (2.22a) 

Clearly, the above interval contains the MP interval and unless s or t equals 0, it is 
strictly larger. In terms of the ratio r, expression (2.22a) takes the form 

g(r) € /[— 3min(l, |r|),3min(l, |r|)], (2.226) 

i.e., the graph of g lies in the shaded region of Fig. 4a. 
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(a) Hyman’s extension (b) M2 region 

Figs. 4. Extensions of the monotonicity region. 


The following limiter function, which was introduced earlier by Van Albada, Van 
Leer, and Roberts [25], satisfies condition (2.22): 

*( r ) = ( 2 - 23 “) 

or 

s 2 t -I- st 2 

G(s,t)= S -^±^~. (2.23 b) 
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Observe that as r tends to oo, g(r ) tends to 1. It can be shown that for r > 0, the above 
limiter lies between the minmod and the harmonic mean limiters. Therefore, it shares 
the property of producing curves that are “too flat”. On the positive side, adjacent to 
discontinuities, this limiter approximates /'(x,) better than any limiter mentioned in 
the previous subsection, as shown in the following example. 

Let x, = i + 1/2 for i = —3, . . . , 2, and 

fi = — 5 — x,- for X{ < 0, fi = 5 — Xi for x* > 0; (2.24) 

i.e., the data lie on the lines y = —5 — x for x < 0, and y = 5 — x for x > 0. Any 
limiter function in the previous subsection gives /_ i = /o = 0. The exact solutions are 
/'(— 1/2) = /'( 1/2) = —1. Using the above limiter, one obtains /_ i = /o = —72/82 
which are much closer to the exact solutions. In the context of conservation laws, 
however, this limiter, similar to the minmod, smears discontinuities. 

The algorithm corresponding to (2.22a) is obtained by replacing fi by the median 
of fi, — 3min(|s|, \t\), and 3min(|s|, |t|). The result is 

fi <- sgn(/,)min(|/,|,3|s|,3|i|). (2.25) 

The above algorithm is very practical due to its simplicity and a slight gain of accuracy, 
see [6], [17], [19]. Theoretically, however, it is only as accurate as the MP algorithm: 
in the example at the beginning of this subsection where the data lie on a parabola, it 
gives the same result, i.e., it still produces only first-order accurate fi near strict local 
extrema. Moreover, it may create oscillations as shown in the next examples. 

Let x,- = i for i = 0, . . . , 4, and /, = i for i — 0, . . . , 3. Let f 4 = k > 4, then the 
data are increasing. The quartic through these points and its slope at x 2 are (see also 
Hyman’s fourth-order finite difference formula in [17]) 

p(x) = x + ^px(x - l)(x - 2)(x -3), A = • ( 2 - 26 ) 

If k = 20, then / 2 = —1, and Hyman’s algorithm (2.25) does not change / 2 . We have a 
strictly negative slope in spite of the fact that the data axe strictly increasing. In fact, 
even the MP algorithm (2.9) gives an “unrealistic” albeit monotone result of / 2 = 0. 

Note that in the above example, if one uses the parabolic formula (2.3), then / 2 = 1, 
and both the MP and the Hyman algorithms work well in this case. In general, the 
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combination of the parabolic formula and Hyman’s algorithm may still create unwanted 
oscillations, as in example (2.24) where the parabolic formula yields /-i = fo = 4, and 
the Hyman’s algorithm then gives /_ i = / 0 = 3. The results are of the wrong sign 
compared to the exact solutions. 

The above examples, especially (2.26), show that one has to be cautious in ex- 
tending the MP interval (2.7) as well as in using higher-order methods to define /,. 
For limiter functions, an extension of the monotonicity constraint for r < 0 can be 
obtained as follows: first, the region must contain the negative r-axis; second, due to 
the symmetric property (2.11), g(- 1) = 0; third, in example (2.24), /_ i and /<, should 
be negative, i.e., for -1 < r < 0, g(r) should be negative, and for -oo < r < - 1 , 
g(r) should be positive. Together with (2.12), these conditions imply that the graph 
of g lies in the region bounded by g = 0, g = 1, g = r, and r = —1, as shown by the 
shaded region for r < 0 in Fig. 4b. For r > 0, the condition is the same as the MP2 
condition in the previous subsection. This enlarged monotonicity region is named the 
M2 region shown by Fig. 4b. Limiters in this region preserve monotonicity in the sense 
that if the data are monotone in [xj,x i+ i], then so is the interpolant. This follows 
because monotonicity of the data in [x,-, ii+i] implies s<_i/ 2 ? 5 i+i/2> 5 i+3/2 8X6 °f 

the same sign. Adjacent to strict local extrema, these limiters do not simply define /, 
as 0; however, they are still only first-order accurate. Two examples of such limiters 
are Van Albada’s limiter (2.23), and 


g(r) = median (1, r, — r — 1) 

= 1 + minmod (r — 1, — r — 2), 


(2.27) 


or 

G(s,t) = median (s,t, —s — t) 

= t + minmod ( s — t,—s — 2 1). 

Note that for r > 0, the above limiter is identical to the minmod. One can also combine 
it for r < 0 with any limiter in subsection 2.2 for r > 0. 

Recently, in [6], a second-order accurate constraint for /, was developed based on 
Hyman’s algorithm. This method essentially detects whether the value /, is at the left 
or the right of a strict local extremum; in each case, a different constraint is enforced. 
Like Hyman’s method, it may produce derivatives of the wrong sign under extreme 
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conditions. In the case of example (2.24), the results are /- 1 — fo — 3. As for example 
(2.26), the result is identical to that of the MP algorithm, which is “unrealistic”: f 2 = 0. 
In the next section, we present a systematic approach in which all data points, except 
those at the boundaries, are constrained by the same conditions which extend that 
of De Boor and Swartz. These extensions produce third-order accurate interpolants, 
allow nonmonotone interpolants for nonmonotone data, and behave well in the extreme 
conditions of the above examples. 


3. Third-order accurate monotone interpolants. To obtain uniform third- 
order accuracy, it is logical to consider the parabola Pi{x ) = (x, p,(x)) through the 
points Fj-i, Fi and F,+i. Denote the second divided difference by d t , i.e., 

ii+Ul ~ tizlll (3.1) 




*^ 1+1 *** 1—1 


Then from the Newton formula, 


Pi+$(x) = fi + Si+i/ 2 (x — %i) + di+o(x — X{)(x x<+i) (3-2) 

where $ = 0 or 6 = 1. If {/,} are third or higher-order samples of a C 3 function /, 
then for all x £ [ij-i, 

Pi(x) = f{x) + 0(h 3 ), p'i(x) = f'(x) + 0(h 2 ). (3.3) 


3.1. Nonoscillatory parabolic interpolation. In the interval [x t , x i+1 ], even 
when the data are monotone, neither pi nor p t+ i are necessarily monotone. However, a 
monotone quadratic Pi +i/ 2 can be obtained by using Harten and Osher s nonoscillatory 
interpolation [15]. We present this method in our geometric framework below. Observe 
that the two parabolas P , , P,+i , and the line segment F, F,+i have the points F, and 
F I+ i in common. Any two of these three curves are either identical or have no other 
point in common. Let P 1+1 / 2 be the curve which lies between the other two, i.e., the 
median of the above three curves, see Figs. 5a, b. Note that the equation for Fi F , 1 is 
given by (3.2) with di+g = 0. Let dj+i /2 be the median of 0, d, , and d»- n» 

d,+j /2 = minmod(di,dj+i), (3.4) 

then the equation for F +1 / 2 is given by (3.2) with 6 = 1/2. 
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Figs. 5. The parabolas -Pj+i/ 2 - 


Before proving the accuracy and monotonicity of p,-+i/ 2 , we review some facts 
for the parabola P(x ) = (x,p(x)) through F,, F t+ i with second-divided difference 
(coefficient of x 2 ) d: 

p(®) = ft *f" •Sj- •£«) + d(x (E,')(x 3-1 + 1 )> 

p'(x) = s i+1 /2 +d(2x-Xi -Xj+i), (3.5a) 

p'(x,0 = -Si+ 1/2 + d(n - x t+ i), p'(xj +1 ) = 5,- +1 /2 + ^(^.+1 — Xi). (3.5 b) 

Denote 

x i+ 1/2 = (®i + *i+l )/2, /i+i/2 = (/. + /t+l)/2i (3.6) 

then the two tangents to the parabola P at F, and F l+ i intersect the vertical line 
x = x i+1 / 2 at the same point B{+ 1 / 2 = (xj +1 / 2 , fei+ 1 / 2 ) where 

bi+ 1/2 = /i+ 1/2 _ -d(Axj + i/ 2 ) 2 . (3.7) 

From (3.5a,b), and (3.7), p is monotone in [x,,x,+i] if and only if any one of the 
following conditions holds: 

Ml < (3.8a) 

AXj +1 /2 
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p'{xi) € J[0,2sj + i/ 2 ], 

(3.86) 

p'(x«+i) € /[ 0 , 2 si + i/ 2 ], 

(3.8c) 

bi+ 1/2 € J[/i>/»+ il- 

(3.8d) 


To show accuracy for p 1+1 / 2 , observe that since d,+i/ 2 is the median of 0, d, and 
d;+i, it lies between di and d,-+i. Equations (3.2) with 8 = 0, 1, and 1/2 imply that 
p i+1 / 2 (x) lies between Pi(x) and p,+ i(x) for all x. Similarly, equations (3.5a) with 
d = di, dj+ 1 / 2 , and d<+i imply p ' i+ 1 / 2 ( x ) i s median of s i+1 / 2 , p((x), and p'i+i(x); 
therefore, p , j+i/ 2 ( a ') lies between p',{x) and p[_|_ j (x ) for all x. As a consequence, it follows 
from (3.3) that for all x € [x,-,x,+i], 

Pi+ 1 / 2 (x) = fix) + 0(h 3 ), p' +1/2 (x) = /'(x) + 0(h 2 ). (3.9) 

We turn now to monotonicity. Suppose the data are monotone in [xj,Xi+i]. It 
will be shown that the same is true for p l+ i/ 2 ; therefore, (3.8a,b,c,d) hold for p,+i/ 2 - 
Consider for the moment the case /, ^ /,•+ 1 . By using a linear change of coordinates 
as in the proof of Lemma 1, we may assume x,' = 0, x,+i = 1, and fi = 0, fi+ 1 = 1. 
The data are thus increasing, and fi - 1 < 0, /,+ 2 > 1. Because s,+i/ 2 = 1, i/ 2 > 0, 

and x I+ i - Xj_i > 1, (3.1) implies d, < 1. Similarly, d i+l > -1. Since d i+1 / 2 is the 
median of 0, di, and d,+i, 

—1 < di + 1/2 < 1 - 

The above expression and (3.5a) imply 


0 <p',+i/ 2 (*) < 2 ( 3 - 10 ) 

for x € [0,1]. Consequently, p, + i/ 2 is strictly increasing in [0,1]. The case /, = f i+ j 
is trivial since the quadratic p,.fi/ 2 reduces to a constant function. This completes the 
proof. 

Additional conclusions can be drawn if the data are monotone in [x,-, Xj.fi]. Indeed, 
we will show that 

P»-l/2(*»)Pi+l/2(*i) > 0, p[ + i/ 2 (x,+l)p[ + 3/ 2 (x,+l) > 0. (3.11) 
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In the above proof, since d, < 1 and d,+i > — 1, equations (3.5b) imply 


> 0, p[+i(zi+i) > 0. 

The facts that 3 i _ 1 / 2 > 0 and s i+3 / 2 > 0 then imply 

Pi- 1/2( X ‘) — Pi+ 3/2( X «+l) — 


Expressions (3.11) follow from the above and (3.10). The case of /, = /*+i is again 
trivial. 

Furthermore, if the mesh is uniform, then in the proof of (3.10), x ,+ \ — Xj_i = 2, 
and the monotonicity of the data in [x,-, Xj+i] implies 


M1+1/2I < 


l^H-1/2 1 
2Ax i+1 /2 * 


(3.12a) 


Pi+1/2( X ) ^ I 


1 3 

2 5 i+l/2) 2 a «+l/2 


(3.126) 


for all x 6 [x,,xj+i]. The equal sign in (3.12a) and the closed interval in (3.12b) are 
necessary for the case /, = /,+i . 

Similar to the definition of s*, let 


ti = minmod [p'_ 1 / 2 (a:.),Pi + i/ 2 ( ;r «)]- (3.13) 

From (3.9), t, approximates /'(x,) to second-order accuracy. Note that the UNO 
scheme [15] is defined by setting /, = t t . Loosely speaking, this scheme favors mono- 
tonicity over accuracy by defining fi to be as close to 0 as possible from the above left 
and right slopes. It can also be considered as a high-order extension of the minmod 
limiter. In the next subsection, we present more accurate schemes which extend all 
limiters. Recall that in the second-order case presented in section 2, limiter functions 
simplify the algorithms by simultaneously satisfying the three conditions of symme- 
try (2.11), accuracy (2.12), and stability (continuity), together with the condition of 
monotonicity (2.13). In the third-order case below, the reverse is true: it is simpler to 
separate the former three conditions from the monotonicity condition. Consequently, 
for third or higher-order methods, limiters only need to satisfy these three conditions. 
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Uniform third-order monotonicity constraints are enforced in a separate step. In the 
next three subsections, we present three different third-order constraints. 


3.2. M3 (monotone third-order) methods. To relax the monotonicity con- 
straint of De Boor and Swartz, (2.6) is replaced by 


Uei 


0,3s i+ 1/ 2, nPt+l/2( X ») 


/i+i e I 


0> 3^1+1/25 fjPi+l/2( X »+l ) 


(3.14) 


The above condition preserves monotonicity, i.e., if the data are monotone in [xi, Xj+i], 
and (3.14) holds, then the interpolant (2.2) is monotone in [x,-, x, +1 ]. Indeed, by ap- 
plying (3.8b, c) to p i+ 1 / 2 , 

2Pi+i/2( x »')> 2^*+ 1 / 2 ^ z * +1 ^ ^ / [0, 3sj+i/2]. (3.15) 


Therefore, the intervals in (3.14) reduce to those in (2.6), and monotonicity follows 
from Lemma 1. Next, expression (2.6') is replaced by 


fie I 


0, 3s;_i/ 2 , «P<-l/2( x ») 


fie I 


0, 3s,+i/2j 2P«+1/2( X *) 


(3.14') 


or fi lies in the intersection of the above two intervals. The algorithm can be simplified 
by observing that the above intersection contains the interval I [0,3s,, |i,] since lies 
in each of the intervals in (3.14'). Consequently, the MP constraint (2.7) is replaced by 
the following condition, which implies (3.14'), 


he I 



(3.16) 


Observe that if the data are monotone in [x^_i , x,] and [x,, x,+i], then the above interval 
is identical to the MP interval. If the data are no longer monotone in [x,_!,Xi] or 
[x,-,x; + i], however, it may be strictly larger. The two end points of the interval in (3.16) 
can be obtained by taking the maximum and the minimum of the three arguments. 
Further simplification is made possible by showing that s, and t, are of the same sign. 
Indeed, P ;_ 1 / 2 ( x «) 1S median of Si_i/ 2 , p'j(xi)> and p(_j(xj) ; p| +1 / 2 ( x «) median 

of s i+ i/ 2 , p',( x i)) and p' i+l (xi)] consequently, 


Pi — l/2( x ») e I — 1 /2 t PtC- 1 -* )]? Pi+1/2( X *) e I [ 5 i+l/2iPi(* l '«)] > (3-17) 
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Because pj(x,) € I [s,-i/ 2 , Si+ 1 / 2 ] by (2.3), the above expressions imply 

m -1/2( X «))P»+1/2( X *)] C I[Si_\j2, 5,+i/ 2 ]. 

Since U = median (0,p-_ 1/2 (x,),p| +1/2 (z,)), and = median (0,«i_ 1/2 ,« i+1/2 ), 

/[0,<i]D/[0,4 

Consequently, s, and t, are of the same sign, and if t, = 0, then s; = 0. Expression 
(3.16), therefore, is identical to 

fie I O.sgn (t,)max(3|s,|,||t 1 |) . (3.16') 

Note that if ti = 0, it does not matter whether sgn(t,) is defined as —1, 0, or 1 in the 
above expression since the interval reduces to the point {0}. 

We summarize the resulting algorithm below. 

Theorem 1. Let s t+l / 2 , di be the first and second divided differences as in (2.1), 
(3.1). Let Si, d i+ 1/2 be the corresponding minmod defined by (2.8) and (3.4). Recall 
the slopes of the parabolas (3.5b) and expression (3.13), 


Pi — 1 / 2 ( x ») = s i- 1/2 + dj_ 1/2 (x, - X,.!), 

Pi+1/2( X ») = 5 i+ 1/2 + ^»+l/2( x « — x i+l)j 
ti = minmod [p{_ 1/2 (*i),Pi +1/2 (*i)]- 

In addition, for m + 2 < i < n — 2, let the original fi be defined by any continuous 
limiter function with property (2.12) applied to the left and right slopes p'_ 1/2 (x,) and 
P'i+i/ 2 ( x «)> e -8-’ arithmetic mean, 

fi = 2 [P«-1/2( X *) d" Pi+ l/2( a '*)]’ (3.18) 

The original fi can also be defined by any second or higher-order formula for the 
derivative. Let the final fi be defined by 


tmax = sgn (ti) max ^3|si|,||<,-|^ , 


(3.19a) 
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/. «- minmod(/i, i max ). 


(3.196) 


Then the algorithm is stable. It preserves monotonicity, i.e., if the data are monotone 
in [x,-,x,+i] for any i, m + 2 < i < n — 3, then so is the interpolant. Furthermore, if 
f £ C 3 , and {/,} are third or higher-order accurate, the final {/,} are second-order; 
consequently, the interpolant (2.2) is third-order. 


The above method is named M3-A (A stands for average) if (3.18) is used. The 
proof is now straightforward. First, stability follows since each of the above equations 
is stable. Next, monotonicity is preserved since (3.19a,b) imply (3.16), which in turn 
implies (3.14') and (3.14). As for accuracy, observe that the original fi defined by 
(3.18) or by a limiter with property (2.12) is second-order accurate. The interval in 
(3.16') contains the value t,, which is also second-order accurate. Any value which 
lies between the original /,• and t, is therefore second-order, in particular, the final /, 
defined by (3.19a,b). 

Before presenting the next method, note the following. By considering the case 
p|_ 1 / 2 (xj) and p' +1 ^ 2 (xj) of the same sign or of opposite sign, the M3-A method (3.18), 
(3.19a,b) can be simplified to: 


fi = sgn (t{) min 


~\p'i-l/2( X i) + P'i+l/2( X i)l 


max 



If the “Superbee” limiter is used, the resulting M3-B method can be expressed as 
fi = sgn (ti) min 

Observe that in the above theorem, if one sets d, = 0 for all i, then the M3-A and 
M3-B methods respectively reduce to the average and “Superbee” limiters. 

The constraint (3.19a,b) has a five point stencil from i — 2 to i -f 2. By using the 
quartic though these five data points (see formula (4.21)), one can obtain a fourth- 
order accurate /,. As shown in example (2.26), such /,• may have the wrong sign. This 
problem can be solved by enforcing the following condition which is a higher-order 
extension of the accuracy condition (2.12): /, lies between p'._ 1 / 2 (x<) and p|. +1 , 2 (a;«), 
or 

fi 4 - median (/,,p'_ 1/2 (xi),p' +1/2 (x,)). (3.21) 


max (|p'_ 1/2 (x,)|,|p' + i /2 (x,)|) ,max ^3|sj|,^|ti|J . (3.20) 
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The final fi is then obtained by (3.19a,b). This method is named M3-quartic. 

If the original /, is defined by a limiter function which is bounded above by 3/2, 
then (3.19a,b) are redundant because the original fi satisfies (3.16) a priori. 

If the data represent a discontinuous function, one can replace (3.18) and (3.19a,b) 
by applying Van Albada’s limiter (2.23) to P-_j/ 2 (x,), P- +1 / 2 ( x »)> and the resulting M3- 
V method is more accurate near discontinuities. (Incidentally, the graph of this limiter 
has a V shape.) It is also highly accurate near the smooth parts of the data because 
the limiter (2.23) satisfies $r'(l) = 1/2. As for monotonicity, suppose that the data 
are monotone in [xj,x, +1 ], then from (3.11), and the fact that this limiter is bounded 
above by g(l + \/2) = (1 + \/2)/2 < 3/2, one concludes that condition (3.16) is satisfied 
at xi and x,+i . Consequently, the M3-V method preserves monotonicity even though 
(3.19a,b) are omitted. Notice that with this method, condition (3.16) may not hold 
if Pf_i/ 2 ( x »)Pi+i/ 2 ( x «) < 0; however, an extension of (3.16), which is similar to the 
extension from the MP2 (Fig. 2c) to the M2 (Fig. 4b) constraint, holds in this case. 
Lastly, as in Theorem 1, accuracy and stability are trivial. 

If the mesh is uniform, and the original /, is defined by a limiter function which is 
bounded above by 2, e.g., the harmonic mean limiter, then it follows from (3.12b) that 
(3.19a,b) are redundant because the original fi satisfies (3.16) a priori. 

The algorithm of the above theorem is an extension of the author’s SONIC schemes 
developed for conservation laws [16]. The results there show that the algorithm im- 
proves accuracy not only near strict local extrema, but also elsewhere (see subsection 
3.5 below). Note that there are some similarities in the techniques of the SONIC 
schemes and (independently) Yang’s recent artificial compression by slope modification 
[28]. The concepts, however, are very different: Yang’s method modifies /, = t,-, the 
slopes of the UNO scheme defined in (3.13), by an adjustable factor obtained from 
numerical experiments; our methods define constraints which can be applied to any 
high-order accurate derivative. 

3.3. MS3 (monotone simple third-order) constraint. To develop the sec- 
ond third-order constraint, we first define the minmod function with more than two 
arguments. Let Z\ , . . . , z\. be real numbers, and 

a = min(xj,. . . ,z*), /? = max(z 1? . . . , z k ). 
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Let minmod(zi,...,zjfc) be the function which returns the smallest argument if all 
arguments are positive, the largest if they are all negative, and 0 otherwise. Then 

minmod (zi , . . . , 2*) = minmod (a, /?). 

The minmod function can also be expressed as 

minmod ,z fc ) = Amin(|zi|, . . . , |z*|) 


where 

A = ^[sgn Oi) + sgn (22)] 7; [sgn (21) + sgn (23)] . . . -[sgn(zi) + sgn (2*)] 

= ^[sgn(zi) + sgn (22)] - [sgn (22) + sgn (23)] . . . -[sgn(2*_i) + sgn(z*)]. 

Observe that in the above expression for the minmod function, if 2, = 0, it does not 
matter whether sgn (z, ) is defined as — 1, 0, or 1. One can also take successive minmods 
to arrive at the same result, e.g., 

minmod (21, 2 2 ,z 3 ) = minmod [minmod (21, z 2 ), 23]. 

Notice that if A: = 2, these definitions are consistent with those in section 2. The case 
of three arguments can be applied to the MP algorithm (2.9). 

The second method to enlarge the MP interval is similar to the first method, except 
ti in (3.16) is replaced by u,: 

/;€/ 0 , 3 ( 3 . 22 ) 

where 

m = minmod [pJ(z.),Pi-i(*i) 5 P;+i( x i)] • (3.23) 

Expression (3.22) can be simplified by showing that s;, u,, and p[(xi) are of the same 
sign. Indeed, since p[(xi) lies in the interval I [s;-i/2i 5 i+i/2L 

j[o, 5| ]c/[o,p;(xi)]. 

This implies s, and pj(x,) are of the same sign, and if p[(x { ) = 0, then s, = 0. Similar 
conclusions hold for u, and p[(x,). Consequently, 

I 0,3s„^u, = 1 0,sgn(p'(xi))max ^3|s;|,||u,|^j . (3.22') 

The corresponding algorithm is summarized below. 
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Corollary 1. Let Sj+ 1 / 2 , d{ be the first and second divided differences as in (2.1), 
(3.1). Let {pj(x,),p^_i(x,),p(_j.i(xi)} be the slopes of the parabolas given by (3.5b). 
Let Si and u, be the corresponding minmod defined by (2.8) and (3.23). Let fi be any 
second, or higher-order approximation of /'(x,), and for m + 2 < i < n — 2, let the final 

i 

fi be defined by 

W = sgn (pj(xj))max ^3|s,|, , 

ft minmod (/,• , < max ). 

Then the algorithm is stable and preserves monotonicity, and the interpolant (2.2) is 
third-order accurate. 


The stability of the above method is trivial. The proof of accuracy is similar to 
that of Theorem 1 with the observation that Ui is second-order accurate. Next, it will 
be shown that 


’ 3 


' 3 

2 

C I 

3s, • , — ti 


(3.24) 


consequently, monotonicity follows. Indeed, since p|_j/ 2 (^i) € J[p' i (xi),p'_ 1 (xj)], and 

P'i+l/2( X i) € 7 b!(*i)>P;+i (*«)]’ 

I \p'i-l/2( X i)^Pi+l/2( X i)] C /bi( x <)»p5-l(*i).p5+l(*i)]; 


therefore, I[0, <,] D I [0, Uj], and (3.24) holds. This completes the proof. 
The above algorithm can be simplified further by defining 


Ui = min(|p'(xi)|, |Pi+i(*i)l) 

instead of (3.23), and the same conclusions still hold. A drawback of this simplification 
is that in example (2.24), if the original /, is defined by p'(x t ), then the results /_i = 
/o = 3 are of the wrong sign. 

Notice that due to (3.24), the method of the above corollary is not as general as 
that of Theorem 1. The next method, however, is more general than the above two. 

3.4. MG3 (monotone general third-order) constraint. Using the four data 
points Fi-i, Fi, F l+1 and Fi + 2 , the monotonicity-preserving parabola P l+J / 2 defined 
in subsection 3.1 is, loosely speaking, “closest” to the straight line segment FiF{+ j. We 
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begin this subsection by defining a monotonicity-preserving parabola P,+i/ 2 which is 
“furthest” from the line segment. Denote F, +1 / 2 = (^.+ 1 / 21 / 1 + 1 / 2 ) where x t+1 / 2 and 
f i+1/2 are defined by (3.6). Let F£ 1/2 = (x j+1/2 , (or F* 1/2 ) be the intersection 

of the line Fj_iFj (or F;+iFj+ 2 ) and the vertical line x = x i+1 / 2 as shown in Fig. 6, 
then 

/i+l/2 = /« + ^i-l/2^i+l/2, fill/2 — /*+l ~ ^ s i+Z/2^i+l/2 ^ 3 ' 2 ^ 

where the superscripts L and R stand for extrapolation from left and right, respectively. 


*Vl/2 / 

c 


<,^ M i+l/2 = irR i+l/2 

1 /2 


: / 


X 


fi-i 

1 1 — 

^i+2 * 

1 i > 


1 1 1— 1 L 

0 12 3 4 

X 

Fig. 6. The parabola P,+ j/ 2 . 


Let P L (or P fl ) be the parabola through F,, F,+i, whose tangent at F; (or F,+i) is the 
line F,_iF, (or F J+ iFi+ 2 ). Similar to the definition of P i+1 / 2 , let P,+i/ 2 be the median 
of the parabolas P L , P R , and the line segment F,F t+1 . Observe that d L , d R , which 
are the second divided differences of P L , P R , can be obtained by expressing the slopes 
of the tangents Fi_iF,, F l+ iFi+ 2 in the form (3.5b), 

JL = ^+ 1/2 - 3 .- 1/2 jr _ + 3/2 ~ Sj+i/2 (3 26 ) 

Ax <+1 / 2 ’ Ax i+1 / 2 


Let 


d( + 1/2 = minmod(d L ,d R ), 


(3.27a) 


then di+i/ 2 is the second divided difference of Pj + i/ 2 , and from (3.26), 


d 


•+1/2 


fi+1/2 

Ax i+1/2 


(3.276) 
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where 

^1+1/2 = minmod(5j +1 / 2 — ■s»-i/2>- s i+3/2 — • s »+i/2)* (3.28) 

Observe that by (3.1) and (3.26), di lies between 0 and d L \ <f,+i lies between 0 and d R ; 
consequently, by (3.4) and (3.27a), 

^i+i/2 £ 7[0,Jj+ 1/2]- (3.29) 

Note that if Pi+ 1/2 has a strict local extremum in (x;,Xj+i), then the same is true for 

Pi+1/2- 

Next, define 

fiti /2 = median (/j + i /2 , Z^i/2, /j+ 1/2 ), (3.30) 

then F™ X f 2 = ( x i+i/2> /i+1/2) i s the intersection of the tangents to P I+1 / 2 at F x and 
Fi+ 1- The slopes of these tangents are given by (3.5b), 

Pi+ 1 /2 C* 1 " * ) = s i+l/2 ^i+l/2> (3.31a) 

Pi+l/ifei+l ) = s i+l/2 4* Si+l/2- (3.316) 

Using definition (3.30), observe that if f+[i/ 2 is a strict local maximum , i.e., f t +. 1 / 2 > fi 
and f^h /2 > /»+i, then from (3.25), the data {/*} have a strict local maximum in 
[X|, i.c., 

/«>/«- 1, and fi+i > fi+ 2 , (3.32) 

see Fig. 6. A similar statement holds for the minimum. Consequently, if f^h/ 2 is a 
strict local extremum , then the data {/;} have a strict local extremum in [xj,x l+ i]. In 
particular, if the data are monotone in [xi,x,-+i], then f^h/ 2 is not a strict local ex- 
tremum, i.e., f^i/2 lies between /, and fi+i- Therefore, by (3.8d), Pi+1/2 is monotone 
in [xj,Xj + i], and by (3.8b,c), 


Pi+ 1 /2 ( *^ * ) ^ 7 [0, 2Sj_)_j/2]> Pi+1 /2(' I, +1 ) ^ 7[0,2s, +1 / 2 ]- 


(3.33) 


Similar to (3.14), (3.14'), we can now relax the constraints (2.6) and (2.6') to 

3 


/, € J 


3 

0 ) 3 s i + 1 / 2 , ^P'i+\/2( x i) 


, and fi+i 6 I 


<yPi+l/2( X i+l) 


I (3.34) 
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fi € I 0,33,_i/ 2, |pi_i/ 2 ( x i) ) fi € I 0,35j +1 / 2 > 2 ^+i/ 2( x «) * (3.34') 

The intersection of the two intervals in (3.34'), denoted by [a*,/?*], is nonempty (con- 
tains 0) and can be obtained by 

3 I 3 I 

af = min 0, 3^j_i/ 2 > 2^*+i/ 2 (* r ’ ) ’ = max 0j3sj_i/ 2 , gPj+iM 1 *) > (3.35a) 

o I 3 ] 

af = min 0, 3sj +1 / 2 > 2 $+i/ 2( x ») > = max 0, 3s,+i/ 2 , -p' +1 / 2 (x,) > (3-356) 

ai = max (of, of), /?, = min(^f ,(3?). (3.35c) 

Notice that a simplification similar to (3.16) does not work in this case because it may 
cause a loss of accuracy. To bring fi into the interval [ft,, /?t ] , again we replace f , by 
median (/j, c*i, /?i ): 

fi *- fi + minmod (a,- - Pi - fi). (3.36) 

Next, we prove that the above algorithm gives a second-order accurate fi if the 
original fi is calculated by a second or higher-order method. In fact, we will show that 

[aiji] D I 0,3 a i+1/2 ,|<i , (3-37) 

which implies ti G [a,-, /?<], and accuracy follows as in the proof of Theorem 1. To show 
(3.37), it suffices to show that the intervals in (3.34) contain the corresponding ones in 
(3.14), i.e, 

I 0,3s i+ i/2,|pi + i/ 2 (®i) 3 7 0,3s i+1 /2,-p' +1 / 2 ( x i) * (3.38a) 

3 nr 3 1 

I 0,3s l+1/2 ,-p' +1/2 (x j+1 ) D I 0, 3s, + i/ 2 , -p' i+1 ^ 2 ( ;r «+i) • (3.386) 

Indeed, from (3.29), 

p!+i/ 2( x «) e 7[s t+ i/ 2 ,p- +1/2 (x,)], P;+i/ 2 ( x «+i) 6 7[5.+i/2,Pi + i/2( x «+i)]’ (3-39) 

which imply (3.38a,b). This completes the proof of accuracy. 

The above method is summarized below. 
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Theorem 2. Let a,, /?,• be defined as in (3.35a,b,c) using (3.28), (3.31a, b), and let the 
original {/, } be calculated by (2.3) or some higher-order method for m + 2 < i < n - 2. 
Let the final { fi } be obtained by (3.36). Then the algorithm is stable and monotonicity- 
preserving in [x,,x :+ i], m + 2 < i < n - 3 . Moreover, if f € C 3 , and {/,} are 
third or higher-order accurate, then the final {fi} are second-order; consequently, the 
interpolant (2.2) is third-order. 

We conclude this subsection with the following remarks: 

Observe that a simplification similar to (3.16) may cause a loss of accuracy since 
in spite of (3.39), the interval l[0,minmod(p'_ 1/r2 (xj),p| +1 ^ 2 (xi))] may not contain 

Due to (3.37), the above method is more general than the previous two. As will 
be shown in the next subsection, the algorithm still provides “plenty of room” for fi 
near strict local extrema if the coefficients 3/2 in expression (3.35a, b) are replaced by 
1. 

Suppose d,+i /2 is an accurate, second divided difference obtained by evaluating, 
e.g., the second derivative at x = Xi+ 1/2 of the cubic through Fi- 1 , F,, F.+i , and F,+ 2 - 
Then, in order that the corresponding parabola A+ 1/2 is monotonicity-preserving, one 
replaces d, +1 / 2 by minmod (d,+i / 2 , d,+i/ 2 ) ■ Note that Theorem 1 with d i+ i/ 2 replaced 
by dj +!/2 is still valid. 

In practice, one can often avoid the calculations in the above theorem by using 
one of the following two tests: 

1. Checking if fi lies in the MP interval by testing if 

Mfi~ 3s.) <0. 

If this condition holds, one moves on to the next point. Otherwise, one goes 

through the calculations of the above theorem. 

2. Checking if the data have a strict local extremum in [x;,x,-+j] by testing if 

(fi ~ /i-l)(/«+l - fi+ 2 ) > 0- 

If the data have neither a strict local extremum in [x,_i,Xj], nor one in [x;,Xj+i], 

we simply use the MP algorithm (2.9). 
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3.5. More on the accuracy of third-order methods. Observe that if the 
data axe on a parabola, then any of the above three methods M3, MS3, and MG3 
reproduces the parabola (if the original {/,} axe calculated by (2.3) for the last two 
methods). However, if the data are on a cubic, these methods may not reproduce the 
cubic even if the original {/i} are exact. We show an example of this case below. 

Let Xi = i for i — —2, . . . , 2, and 

f—i = f 0 = /i = 0, f—2 = -6, and / 2 = 6, (3.40) 

i.e., the data are on the cubic f(x) = x 3 - x. Clearly, the parabola P 0 through F_i, 
F 0 , and F x is identical to the x-axis. Since the interval [ao,/?o] reduces to {0}, any 
of the above three algorithms yields /o = 0. The exact solution is f (0) — —1, and 
the algorithms do not reproduce the cubic. If h is small enough, however, the MG3 
constraint reproduces the above cubic as shown below. 

More generally, suppose the mesh is uniform, the data axe fourth or higher-order 
samples of a C A function /, and the original {/,} are third or higher-order accurate, 
then when do the above third-order constraints cause a loss of fourth-order accuracy? 
As shown at the end of subsection 2.1, at the part where f'(x) ^ 0, the MP algorithm, 
and therefore all third-order algorithms, does not alter provided h is small enough; 
consequently, accuracy is preserved. Next, we show that adjacent to a strict local 
extremum at x = x* where /"(**) / 0, the M3 and MS3 constraints are essentially 
identical, and they may yield only second-order accurate /, ■ The MG3 constraint, on 
the other hand, does not alter fi in this case, and thus, preserves accuracy. Indeed, 
for any xi, after a linear change of coordinates, we may assume x t = 0, and the Taylor 
series expansion of the exact solution at x = X{ is of the form 

/(x) = sx + dx 2 + ex 3 + 0(h 4 ). (3.41) 

Higher-order terms will be omitted below when there is no ambiguity. By evaluating 
/i— 2 ,---,/i+ 2 > ^ using (2.1), (3.1), and (3.5b), 

5j _ 3 /2 =s -3dh + 7eh 2 , s, +3/2 =s + 3dh + 7eh 2 , (3.42a) 

s,_i/ 2 = s — dh + e/i 2 , Si+ 1/2 = s + dh + eh 2 , (3.426) 
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di-i=d — 3eh, di = d + 0{h 2 ), d,+i = d + 3eh, (3.43) 

Pi_i(x.) = s - 2eh 2 , Pi(xi) = 3 + eh 2 , p' i+1 (xi) = s - 2eh 2 . (3.44) 

Since / 0, in a small enough neighborhood of x*, we may assume d 0. The 

following equations (3.45-50) only use the fact that d / 0. Suppose h is small enough 


so that 

|d| > 3|e| h. (3.45) 

If de > 0, then from (3.4), (3.5b), (3.27b), and (3.31a, b), 

di— 1/2 = d — 3eh, ^ 1 + 1/2 = (3.46) 

Pi-i/tM = 3 -2eh 2 , Pi+r/aC**) = 5 + (3-47) 

di— 1/2 = 2d — 6e/i, ^i+ 1/2 ~ 2d. (3.48) 

P«'-i/ 2 ( x «) = 5 + d/» - 5e/i 2 , p' i+1/2 (xi) = s-dh + eh 2 . (3.49) 

Note that since the mesh is uniform, d,±i/ 2 = 2di±x/ 2 . If de < 0, then 

di- 1/2 = d, d, + i /2 = d + 3eh, (3.46') 

Pi — 1 / 2 * ) = ^ 4” ^ > Pi-f 1/2 (■*'•) = ^ 2eh , (3.47 ) 

di- 1/2 ~ 2d, ^i+ 1/2 = 2d + 6eh. (3.48') 

Pi_i/ 2 ( x «) = 5 + d/i + eh 2 , P'i+\ /2 ( x i) = 5 - dh - beh 2 . (3.49') 


FYom (3.44), (3.47), and (3.47'), the M3 and the MS3 constraints are identical up to 
0(h 3 ): 

ti = u, = minmod(s + eh 2 ,s — 2eh 2 ). (3.50) 

They may yield only second-order accurate fi adjacent to a strict local extremum 
because of the following reasoning: if e ^ 0, and s € I [— eh 2 , 2eh 2 ] such that s is 
strictly 0(h 2 ), e.g., s = eh 2 , then ti = it, = 0 by (3.50). By (3.45), \d\h > 3\e\h 2 . 
Equations (3.42b) then imply s; = 0. Consequently, the final /, is equal to 0 by either 
of these two methods, and the error is 0(h 2 ). Note that if |s| > 6|e| h 2 , then the M3 
and MS3 constraints preserve high-order accuracy since they do not alter the original 
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fi due to (3.50). A loss of fourth-order accuracy, therefore, may occur only in a small 
region around a strict local extremum whose length is 0(h 2 ). 

On the contrary, the MG3 constraint preserves accuracy: if h is small enough, then 
up to 0(h 2 ), the interval [»*, /?,-] of the MG3 constraint defined in (3.35a, b,c) contains 
I[s - dh,s + dh] due to (3.42b), (3.49), and (3.49'). Since d / 0, this constraint does 
not alter the original fi = s + 0(h 3 ). The MG3 constraint (and therefore, all third- 
order constraints), however, can cause a loss of accuracy near x where f'(x) = 0, and 
f"(x ) = 0, as shown in the following example, see also [7]. 

Let the data be on the cubic f(x) = x 3 , and let the mesh be uniform with x,_! = 
— c, x,- = 2e where e > 0. It follows that 

/(-I e 3 , /, = 8e 3 , /.+, =125e 3 . 

Since the data are increasing, the interval [a,,/?,] reduces to [0,3s;] = [0, 9e 2 ]. The 
exact solution, f'(xi) = 12c 2 , lies outside this interval, and the final /, has error 
0(e 2 ) = 0(h 2 ). 

We conclude this section by comparing the accuracy of /; defined by different 
third-order methods. Observe that the (lowest-order) error of the UNO scheme (the 
h 2 term of |<, - /'(x,)|)lies between \eh 2 \ and \2eh 2 \ by (3.50). The error of the M3-B 
method also shares this property due to its definition, (3.47), and (3.47'). The parabolic 
formula defined by fi = p'-(xj) has error |e/i 2 |; consequently, it is more accurate than 
the UNO and the M3-B methods. Away from strict local extrema, the error of the 
M3-A method is |e/i 2 /2|, which is roughly less than half that of the UNO, the M3-B, 
or the parabolic methods. Adjacent to strict local extrema, the error is at most |2e/i 2 |. 
Notice that one cannot eliminate the second-order error term by a limiter function due 
to (3.47) and (3.47'). However, this term can be eliminated by using (3.44): 

fi - £ [4p!(*i) +Pi-i(xi) + P;+i(*i)] • 

In fact, if the mesh is uniform, the above formula is identical to the quaxtic formula; 
consequently, the error is 0(h 4 ). In this case, the M3-quartic method consists of the 
above formula followed by (3.21) and then (3.19a,b). It yields a fourth-order accurate fi 
except at strict local extrema and inflection points, where it may degenerate to second- 
order. At inflection points, the loss of accuracy is caused by (3.21), and at strict local 
extrema, by (3.19a,b). 
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In the next section, we extend the monotonicity intervals further and derive fourth- 
order accurate algorithms which reproduce a cubic when the data lie on one. 

4. Fourth-order accurate monotone interpolants. There are three main 
problems in developing fourth-order monotone interpolants which reproduce a cubic 
when the data are on one. First, example (3.40) shows that there can be monotone 
data which lie on a nonmonotone cubic. If the interpolant reproduces the cubic, it 
cannot preserve monotonicity as defined in section 2. Therefore, one has to modify the 
definition of monotonicity. Next, second-order accurate /< can be obtained by using 
formula (2.3), which has a three point stencil consisting of i - 1, i, and i + 1. For 
third-order accurate one has to add one more point to the above stencil, and to 
preserve symmetry, the stencil is enlarged to five points from i 2 to i -f- 2. If f , is 
calculated using the quartic through these five data points, it can be of the wrong sign 
as shown in example (2.26). The second problem is thus to define the original /, to 
third-order accuracy. The final problem is to enlarge the monotonicity interval. 


4.1. Nonoscillatory cubic interpolation. To solve the above problems, we 
introduce the cubics £,'+ i/2( x ) defined for x, < x < x,+i whose graphs {<?,+ i/2( x ) = 
(x,g; + i/ 2 (z))} stay “close” to the straight line segments FiFi+i so that they do not 
create unwanted oscillations. These cubics can be considered as fourth-order extensions 
of the nonoscillatory parabolic interpolant Pi+ 1/2 defined in subsection 3.1. Concep- 
tually, our cubics are somewhat similar to those of the ENO schemes [14] which were 
developed for conservation laws. The key difference is that the ENO cubic interpolation 
is unstable due to a switch which decides whether the stencil is enlarged to the left or 
the right. Our methods are stable because they employ the median function as shown 


below. 

Let q i+ x/ 2 (x) be the cubic defined by the four points F,_i, F,, Fj+i, F t+2 , and let 
e l+ 1/2 be the third-divided difference of {/,}, 


e t+l/2 ~ 


dj-t-i — dj 
Xi_j»2 — 1 


(4.1) 


From the Newton formula, 


9«'-l/2( x ) = Pi( x ) + e i-l/2( x ~ x i-l)( x ~ x i)( x x «+l)) (4.2a) 
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9i+l/2( x ) = Pi( x ) + Gi+\fi(x ~ x i-l)( x x i)( x x »+l)? (4.26) 

where p, is defined in (3.2). Rewrite (4.2a) with index i replaced by i + 1, 

9 i+l/ 2 { x ) ~ Pi+l( x ) + e i+l/l( X — x i)( x — a: i+l)( x “ x «+ 2 )- (4.2c) 

If {/,} axe fourth or higher-order samples of a C 4 function / , then for all x € [xi-i , Xi+2]» 

ft+i/a(*) = /(*) + 0{h% g' +1/2 (x) = /'(*) + 0(h 3 ). (4.3) 

In the interval [xi,Xi+i], the cubic Oi+ 1/2 Soing through F,, F t + 1 is defined by 
the derivatives q' i+1/2 ( x i) and g| +1/2 (a?i+ 1)- We present two different methods to define 

Qi+ 1/2- 

For the first method, in the interval [x t -_i,Xj + i], consider the three curves <?,- 1 / 2 , 
Pi and Q,+i/ 2 , which go through the points F,-\, Fi and Any two of these curves 

are either identical or have no other point in common. Let Qi be the curve which lies 


between the other two, and let 

e, = minmod(ei_i/ 2 ,ej +1 / 2 )- ( 4 - 4 ) 

From (4.2a,b), and (4.4), the equation for Qi can be expressed as 

?i(x) = Pi(x ) + ti{x - X j_ j )(x - X;)(X - Xj + j). (4.5) 

Upon differentiating the above equation, 

q'i( x i-l) = p'i( x i-i) + ei(xi-i -Xi)(x,_i -x l+ i), (4.6a) 

?{(*»•) = p'i( x i ) + c«(x,- - Xi-i)(xi - x i+ i), (4.66) 

q'i( x i+ 1 ) = Pi( x i+i) + e i( x i + 1 “ x,_i)(x i+1 - x,). (4.6c) 


Observe that in the interval [xj,Xj+i], among the two cubics Qi, Qi+i and the straight 
line segment F t 1 , there may not be a median curve as in the case of the parabolas 
Pi, Fi + 1 and FiFi+i shown in the previous section. However, one can take the median 
of the derivatives to define q i+ !/ 2 , 

= median [s i+1 / 2 ,9|(xi),9- +1 (x t )]. (4.7) 
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Using (3.5b) and (4.6), equation (4.7) and (4.7) with a replaced by x i+ i imply 
q' i+ 1 / 2 ( x i) = 5 i+i / 2 — Axi+i/2 minmod [dj + e,(x,- — x,_i ), dj+i + e,+i (x, — Xj+ 2 )] » (4.8a) 
q' i+ 1 / 2 ( x i+i) = 5,+i/ 2 + Ax i+1 /2 minmod [dj + e,(x, +1 -Xj_i),d i+ i + ei+i(x i+ i -x i+2 )j. 

(4.86) 

Remark that if d; and d ,-+ 1 are of the same sign, then the cubic 0*+i /2 is the 
median of the three curves Q > , Qi+i and FiFi+ To show this, observe first that in 
the interval [xj,x,+i], the cubic Q i+ j/ 2 lies between the two parabolas P, and p+i for 
arbitrary d,, dj +1 . Indeed, by differentiating (4.2b,c), 

/2 C*®' * ) = Pi( x i) e i+l/2( x i ~ x i-l)( x i ~ x i+l)> 

(4.y) 

= Pi+ l( x i) + e «+l/2(*« - *»'+l)( x « - X, +2 ). 

Since (x* — x,_i)(xj — x,+ 2 ) < 0, equations (4.9) imply 

[?«'+i/2( x ») — P , «( x i)][9»+i/2( x «) — P«+i( x *)] — 

or 

9i+l/2( X *) ^ I [P«( X «)> Pt'+l ( X «)]' (4-10) 

A similar statement with x, replaced by x; +1 holds. Because any two of the three 
curves Qi + j/ 2 , Pi , and Pi + j are either identical or they do not intersect in (x;, x,+i ), 
Qi+i /2 lies between P, and P,+i, see Fig. 7. Next, since Qi is the median of Q,_ i/ 2 , 
Q.+ 1/2 and -Pn it H es between Q 1 + i/ 2 and P,. Similarly, Q l+ i lies between Qi+ 1/2 and 
P i+1 . Thus if 0 <di < d,+i, or 0 > di > d,+i as in Fig. 7, the parabola Pj is closer to 
P.P.--H than Pj+i and Qi+ 1/2 is identical to Q,. A similar conclusion holds if d l+ i lies 
between 0 and d;. If d, and d i+ i are of opposite sign, however, the cubic Oi+ 1/2 can 
be different from any of the six curves Q,_ i/ 2 , Qi+i/ 2 > Q.+ 3 / 2 , -Pi. P«+i> and P«Pi+i 
which are used to define it. 

The above proof also shows that definition (4.7) is equivalent to 

Qi+i/ 2 ( x i) = median [p| +1/2 (xi),g,-(x,), 9 j +1 (xi)] 

where s i+1 / 2 in (4.7) is replaced by p' i+ 1 / 2 ( x i) defined in (3.5b). Our next definition 
resembles (3.13) and will be used to extend the monotonicity interval, 

f, = minmod [g'_ 1/2 (x,),# +1/2 (x;)]. (4-H) 
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Fig. 7. Qj+i/2 lies between Pi and Pi ti- 

For the second method, Qi+ 1/2 is again defined by the two slopes ?; +1 / 2 (zi)> 
#+i/2( x H-i) where 

u min = m ' n [?j — i/2( ;r *)> 9 l +3/2( 3 '*)] ’ 

u max — — 1 /2 C ^ * ) ’ 9*4-1 /2 (**'*)■> 9f-§-3/2( a '* )] ’ 

( 4 . 12 ) 

1 / 2 (•*"*) = median (s,+i / 2 , u m in , «mw)i 

= Sj^. 1/2 + minmod (l/rnin s i+l/2i u max • s i+l/2)- 

Similarly, <^ +1 / 2 (ar,+i) is defined by replacing x t by x,+i in the above equations. Ob- 
serve that 

I Si+lfa*)] C I [9i-l/2( a '*)’ 9i+l/2( a '*)> 9i+3/2( a '*)]5 

consequently, from definitions (4.7) and (4.12), 

I [ ,s i+l/2t9»-|-l/2('®*)] ^ [ s i+l/2? /2 (**'*)]'» (4.13) 

i.e., q' l+1/2 (xi) and q' i+1 / 2 ( x i) are on the same side relative to s; + i/ 2 and g( +1 / 2 (x,) is 
closer to s i+1 / 2 than 9,- + i/ 2 (^i)- Finally, as in the first method, 

t, = minmod [g'_ 1 / 2 (x 1 ),^ +1/2 (x i )]. (4-14) 

Notice that if the data are fourth-order accurate, then so are {Oi-t- 1 / 2 } defined by 
either of the above two methods. Because these cubics stay “close” to the straight line 
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segments F,F,+i, they do not create unwanted oscillations. However, the derivatives 
of these cubics are discontinuous at {xj}; therefore, they are only C°. We use these 
cubics to define accurate, monotone C 1 interpolants below. 


4.2. M4 (monotone fourth-order) methods. First, we modify the defini- 
tion of monotonicity. The data axe q-monoione in [x,, x,+i] if they are monotone in 
[xi,x l+ i], and 

9i'+i/ 2( x ») ^ -f [0, 3 s, + j/2]j 9»+]/2( ;e *+i ) € I [0, 3s,+i/2] (415) 

where $ +1/2 (xj), 9 ( + i/ 2 (x»+ l) are defined by (4.8). Using Lemma 1, (4.15) implies that 
Q i+ 1 / 2 is monotone in [xj,x,+i]. Observe that if the data are q-monotone, then (4.15) 
with 9 j + i/ 2 ( x i)) 9 !+i/ 2 (‘ e, + 1 ) re pl ace( i by Qi+i/ 2 ( a: *+i) * s a ^ so satisfied due 

to (4.13). 

Next, the original /, is defined by applying any continuous limiter function with 
property (2.12) to the left and right slopes 9 i_i/ 2 ( x «) a 11 ^ i/2(‘ r »)» e -S-> the arithmetic 

mean as in (3.18), 

/« = \ 9,'-i/ 2 ( x «) + 9i+i/2( x <)] • ( 4 - 16 ) 

The corresponding method is named M4-A. If the data represent a discontinuous func- 
tion, the Van Albada’s limiter (2.23), which results in the M4-V method, is preferred. 
Similar to the M3-V method, the monotonicity algorithm (4.20) below is omitted for 
the M4-V method. 

Finally, we relax the monotonicity constraint (3.14) to 

3 ] 

/i € / 0, 3.Sj_(.i/2> ' (4.17 a) 


• r 3 

fi + 1 £ I 0)3s,-|.i/ 2) 2P«+l/2( X ‘ + 1 )’9<-H/2(* C »+l) 


(4.176) 


More generally, we can extend any constraint of section 3. For simplicity, the following 
constraint extends (3.16), 

/, e/ 0,3«„ !(„(,] (4.18) 


where t, are defined by (3.13) and (4.11). 

The fourth-order methods are summarized below. 
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Theorem 3. Let s i+1 / 2 , d», e i+1 / 2 be respectively the first, second, and third di- 
vided differences; let s i} d l+1/2 , e, be the corresponding minmod as in (2.8), (3.4), and 
(4.4). Let {p'_ 1/2 Oi), P' i+ i/ 2 ( x i)} 311(1 U be given by (3.5b) and (3.13); let q- +1 / 2 ( x i), 
q'. +1 f 2 (xi+i), and U be defined by either (4.8a, b) and (4.11), or (4.12) and (4.14). In 
addition, let fi be defined by a continuous limiter function with property (2.12), e.g., 
the average (4.16), or by any third or higher-order formula for the derivative. Let the 

final fi be defined by 

tmin = min (o, 3sj, |t,, , (4.19a) 

Wx = max (0, 3si, | U, t^j , (4.196) 

minmod (i m in — /«>< max — /»)• (4.20) 

Then the algorithm is stable. Moreover, if the data are q-monotone in [x, , - T i+i ] foi 
any i,m + Z<i<n-4, then the interpolant (2.2) is monotone in [xj,x,+i]. As for 
accuracy, if f € C 4 , and {/,} are fourth or higher-order accurate, the final {/,} are 
third-order ; consequently, the interpolant (2.2) is fourth-order. 

The proof of stability is again straightforward. To prove accuracy, observe that 

the second half of (4.3) implies gj-_ 3/2 , <l',- i/ 2 > ?;+ 1 / 2 ’ and ^+ 3/2 at x = x * a PP rox ‘ 
imate the true derivative /'(*<) to 0(6 3 ). Therefore, the same is true for q^, q\, 
q' i+1 because their definitions use the median function. It follows that ?(•_ i/ 2 , ?|+i/ 2 > 
and then t, , are third-order accurate. Since the original f, is defined a limiter func- 
tion with property (2.12), it lies between q\_i( 2 (.Xi) and ?|- +1 / 2 (*i)> and tlius 1S also 
third-order. Because the final fi by (4.20) lies between the original /, and £,, the 
accuracy claim follows. Next, suppose the data are q-monotone in [xi, ^i+i], then 
q' i+1/2 (xi) 6 I[0,Zs i+1/2 }. Definition (4.11) implies t, <E l[0,q' i+1/2 (xi )]. Consequently, 
U £ I [0, 3s, + 1 / 2 ]- Since q-monotonicity implies that the data axe monotone in [x,-, x;+i], 
it follows that | ti G I[0,Zs i+1/2 }. The interval J[0, 3s„ \ uji) is thus contained in 
J[0,3sj + i/ 2 ]- Expression (4.18) then implies fi € I [0,3s; + i/ 2 ]. Similar arguments hold 
for fi+i, and the monotonicity of the interpolant (2.2) follows from Lemma 1. This 
completes the proof. 

Note that instead of the average (4.16), one can use the following fourth-order 
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quartic formula to define the original /<: 

i ,, v . , X, x C«-l/2(*»+a “*«) + e i+l/2(**“*i-2) , A ni'k 

ft = p&i) + (*,- - - *<+i) — — - z 7 ~, : • ' 4 - 21 ) 

x i+2 ^i-2 


As in (3.21), to avoid derivatives of the wrong sign, before applying (4.19a,b) and 
(4.20), one replaces the original /, by 


/, <- median (/j,g-_ 1/2 (ii)7 9|+i/2( x '))- 


(4.22) 


The resulting method is called the M4-quartic method. 

To obtain MS4 constraint, which extends the MS3 constraint, one simply replaces 
<i, U in (4.19a,b) by Ui, it, where U{ is defined in (3.23), and 

Hi = minmod [g'_ 3 / 2 (a;i)>? , i-i/2( x «).?<+i/2( x i)>9!+3/2( :r «)] • 

This constraint, however, does not provide the original /,. If the quartic formula (4.21) 
is used for /;, then (4.22) is necessary to avoid derivatives of the wrong sign. Since 
q'._ i/ 2 ( x ») an< * 9 !+i/ 2( x *) are nee< ied, one might as well use the algorithm of Theorem 3. 

5. Boundary conditions. The above methods apply only to the interior points. 
Near the boundaries, one must use one-sided algorithms. We carry out the arguments 
for the left boundary, and only occasionally, formulas for the right boundary are pre- 
sented to show the symmetry between the two boundaries. For third-order methods, 
f m is defined by the slope of the parabola Pm-t-i at x m , 

fm — Pm+ l( x m) = s m+l/2 "I" < ^m+l( a; m ^m+l), (h-1) 

fn = Pn-l(Zn) = Sn-l/2 + d n -i{x n - X n _j). 

For the quartic or the fourth-order methods, the obvious choices for f m and / m+1 are 
the slopes of the cubic Q m+3 / 2 . It follows from differentiating (4.2b) that 

fm = Sm + l/2 + d m 4 -i(x m — X m + x) + e m + 3 / 2 ( x m — X m + i)(^m — x m+ 2 ), 

f m + 1 = s „ ,-(-1/2 + d m + i(l m + l — x m) + ^m+3/2( x m+l ~ x m)( x m- |-1 ~ x m+ 2)- 

As for the constraints, the second, third, and fourth-order ones presented above 
have respectively a three, five, and seven-point stencil from i — k to i + k, k = 1, 2 or 3. 
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Near the boundaries, the high-order constraints are replaced by appropriate lower-order 
ones, e.g., at m + 2, one of the third-order constraints, or the MP constraint, and at 
m + 1, the MP constraint. Since the definition of monotonicity of the data in [x,, x,+i] 
involves the values fi—i and fi+2i at the left boundary, we cannot have monotonicity of 
the data in [x m ,x m +i]. In this interval, the interpolant is monotone provided condition 
(2.6) is enforced: 

f m <— minmod(/ m ,3s m +i/2), /« «— minmod(/ n ,3s n _i/ 2 )- (5.2) 

If the data are smooth and have a strict local extremum next to the boundary, then the 
above constraints may “clip” the extremum as shown in subsection 2.3. Consequently, 
whenever possible, the boundaries should be chosen at the smooth part of the data and 
away from strict local extrema. 

For the M3 methods, boundary conditions which provide smoother transitions 
with the interior points than (5.1) can be obtained as follows. The task is to define the 
parabola P m . The simplest way is to let P m be identical to P m + 1 , 

d m = dm+\ • ( 5 - 3 ) 

P m can also be defined by the two points F m , F m+ i, and the slope q' m+3 / 2 ( x m)- 
(4.2b), 

Qm+ 3 / 2 { x ) ~ Pm + \ (®) 4" |-3/2 — x m){ x ~ ^m+lX® — x m+ 2 )- 

By expanding q' m+3 / 2 ( x m), and using (3.5b), the above definition for P m is equivalent 
to 

d m = d m+ i + e m+3 / 2 ( x m — x m+ 2 )- (5-4) 

dn = d n _ 1 + e n -3/ 2 ( x n ~ x n- 2 )- 

The third method to define P m is obtained by setting p" m equal to q'm+3/2( Xm ')' or 
equivalently, 

d m = d m+ \ + e m+3 / 2 (2x m — x m .^2 — ^m+i)- (5-5) 

Observe that if the mesh is uniform, (5.4) and (5.5) respectively reduce to 

2 

dm dm-\~ 1 4* 7^ (, dm 1 d,,, 2 ) i 
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dm — d m + 1 + (t/ m +i d m +2')' 


Expressions (5.3), (5.4) and (5.5) for d m correspond to extrapolation to the boundary 
from the interior, with (5.4) and (5.5) being more accurate them (5.3). It is well-known 
that higher-order extrapolation may cause problems. In this case, however, the essential 
quantity is d m + i/ 2 , which is the median of d m , d m +i, and 0. Therefore, d m takes effect 
only when it lies between 0 and d m+1 , and consequently, the above methods behave 
at least as well as the commonly used one-sided interpolation (5.1). To define / m , we 
simply set 

fm — Pm+l/ii 1 ™)' (5-6) 

Clearly, if one uses (5.3), the resulting f m is the same as (5.1). If (5.4) is being used, 
then f m by (5.6) is identical to 


f m = median 


»+l/2 > Pm+1 ( Xm )’ ^m+3/2( Xm )j • 


(5.7) 


Note that to assure the monotonicity of the interpolant in [x m , *m+i] and [i m +i,i m + 2 ], 
the MP constraint must be enforced at m + 1, and (5.2) at m. 

Finally, for the M4 methods, the first task is to define Q m +\/2- Let 


7 = /[x m , 


•• jZ’m-f <l] 


^m+ 5/2 — e m + 3/2 
4 


(5.8) 


and let r m + 2 (z) be the equation of the quartic interpolating f m , ■ ■ ■ , fm+ 4- Then by 

(^) 

applying the Newton formula, and setting e m+1 /2 = r m+2( X m+l/2)/3!, 


Cm+1/2 — e m+3/2 4" 7 (^m d" (Em+l £m+2 (Tm+3)- 


(5.9) 


Next, we define Q m +i/2> 

q'm+l/ 2 ( x rn) = median [s m + 1 / 2 , ?m + i/ 2 (*m), <7m+3/2( x m)], (5.10) 

and similarly, 9j n+ i/ 2 ( ;r m+i) is defined by replacing x m by x m+1 in the above equation. 
Consequently, at i — m-f-1, and i — m- (-2, there is enough information for the algorithm 
in Theorem 3. At the left boundary, set 


fm — 9m + l/2( Xrn )‘ 


(5.H) 
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fm = 
fm + 1 


( 6 . 16 ) 

(6.1c) 


0. Numerical examples. The methods presented above can be separated into 
two types: methods which combine /, defined by an interpolation formula with the 
MP, M3, MS3, MG3, or M4 constraint belong to type I; those which yield the slopes 
on the left and the right side of i before employing a limiter function belong to type 
II. For type I methods, besides the parabolic (2.3) and the quartic formula (4.21), we 
also test Hyman’s fourth-order finite-difference formula [17], 

i _ ~/»+2 + S/,+i — 8/i-i + / «— 2 (6.1a) 

— x,-|-2 + 8xi+i — 8x;_i + X {— 2 

If the mesh is uniform, the above formula is identical to the quartic formula. For 
convenience, we refer to (6.1a) as the FD4 formula. At the left boundary, the following 
cubic formulas are used: 

— 22 fm + 36/m+l — 18/ m +2 + 4/ m +3 
-22x m + 36x m+ i - 18x m+2 + 4x m+3 ’ 

-2 fm - 3/m+l + 6/m+ 2 ~ 2/ m +3 
2x m 3x m+ i -|- 6x m -|_2 2x m .(-3 

Formulas for the right boundary are obtained by replacing m+ by n— in the above 
equations. Type I methods to be tested are the combinations of the parabolic, the 
quartic, or the FD4 formulas with one of the constraints MP, MS3, MG3, M3, or 
M4. Type II includes methods which employ limiter functions. The eight limiters 
to be tested axe: the minmod (2.14), abbreviated by MD; the Van Albada (2.23), V; 
the average (2.19), A; the “Superbee” (2.17), B; and occasionally, the harmonic mean 
(2.15), HM; the Fritsch-Butland (2.16), FB; the limiter (2.20), AR (average rational); 
and the limiter (2.21), AQ (average cubic). Combining these limiters with the M3, M4, 
or no constraint yields 24 methods. For the M3 methods, we use the boundary condition 
(5.5), and for the M4, (5.10) and (5.11). Since there are 39 methods altogether, for 
each data set below, we can only present some selected results for these methods. 

The data sets include both rough and smooth ones. The rough data sets, which 
are monotone, are used to test the “visually pleasing” or “too flat” properties of the 
methods. The smooth data sets are nonmonotone. They provide information on accu- 
racy. 

6.1. Monotone data. The first data set appearing below is an example from 
Akima, [1], [9], [10], and the second is the Fritsch-Carlson RPN 14 data [10], [17]. 
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Akima data 

Fritsch- Carlson data 

X 

f 

X 

f 

3. 

10. 

7.99 

0 . 

5. 

10. 

8.09 

2.76429E-5 

6. 

10. 

8.19 

4.37498E-2 

8. 

10. 

8.7 

0.169183 

9. 

10.5 

9.2 

0.469428 

11. 

15. 

10. 

0.943740 

12. 

50. 

12. 

0.998636 

14. 

60. 

15. 

0.999919 

15. 

85. 

20. 

0.999994 


Figs. 8, 9 and 10, 11 show the respective interpolation curves for the Akima and 
the Fritsch-Carlson data. Since these data Eire monotone, the M3, MS3, and MG3 
constraints axe identical to the MP constraint. 

The methods presented in Figs. 8 and 10 are: (a) parabolic, (b) MG3-parabolic, 
(c) FD4, (d) M4-quartic, (e) MG3-FD4, (f) M4-A, (g) harmonic mean, and (h) Fritsch- 
Butland. Note that the interpolants on the right are generally more “visually pleasing” 
than those on the left. The interpolants in Figs. 8(a,c) of the parabolic and FD4 
methods are oscillatory. Since the mesh is irregular, /„ at the right boundary of the 
FD4 method has the wrong sign in Fig. 8(c). Although the monotonicity constraint 
corrected f n to 0 in Fig. 8(e), it is still unacceptable. Consequently, the FD4 method 
(6.1) should only be used when the mesh is smooth. Fig. 8(g) shows that the harmonic 
mean limiter, as observed by Fritsch and Butland [9], produces interpolants which 
are “too flat”. Four of the methods which yield “visually pleasing” interpolants are 
presented in the right column: the constrained parabolic, constrained quartic, M4- 
A, and Fritsch-Butland. Except for those of the boundary, similar remarks hold for 
Figs. 10. 

In Figs. 9 and 11, results for well-known (M2) limiters are presented on the 
left, and the corresponding M3 methods are on the right: (a) minmod, (b) UNO, (c) 
Van Albada, (d) M3-V, (e) average, (f) M3-A, (g) “Superbee”, and (h) M3-B. The 
interpolants of the minmod, Van Albada, UNO, and M3-V methods, shown by Figs. 9 
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and 11 (a,b,c,d), are “too flat”. Although the M3 methods are more accurate, the 
results are essentially the same as those of the M2 methods since the data axe rough. 
The interpolants for the AR and AQ limiters are not shown because they are essentially 
the same as those of the average. Methods which yield “visually pleasing ” interpolants 
are presented at the bottom half of each figure: the average limiter, “Superbee” , M3- A, 
and M3-B. 

For rough and monotone data sets, “visually pleasing” interpolants axe produced 
by the following methods: the constrained (MP, MS3, or MG3) parabolic; the average, 
Fritsch-Butland, and “Superbee” limiters; the M3-A, M3-B, and M3-quartic; and the 
M4-A, M4-B, and M4-quartic. 

6.2. Nonmonotone data. Following Hyman [17], we interpolate the function 

_ e -x 7 on a uniform mesh. Note that the data are smooth and have a local 
maximum. The results for N = n - m = 8 mesh intervals are presented in Figs. 12- 
14. Since the mesh is uniform, the FD4 formula gives identical results as the quaxtic 
formula, and the MP-parabolic method is the same as the average limiter method. 

Fig. 12 shows the results for the domain [-2.8, 3.6]. The methods presented are 
of type I: (a) parabolic, (b) MG 3-parabolic, (c) MP-parabolic, (d) MS3-parabohc, (e) 
quartic, (f) MG3-quartic, (g) MP-quaxtic, and (h) MS3-quaxtic. Figs. 12(c,g) show that 
the MP constraint clips the strict local extrema, while the MS3 and MG3 constraints 
have no effect in this case. Observe that in Figs. 12(a,e), near the boundaries, the 
parabolic and the quartic methods yield nonmonotone interpolants, while all the con- 
strained methods preserve monotonicity. Accurate monotone interpolants are produced 
by the methods presented in the right column: (b) MG3-parabolic, (d) MS3-parabolic, 
(f) MG3-quaxtic, and (h) MS3-quartic. 

Figs. 13 are identical to Figs. 12 except the domain is shifted to the left: [-2.9, 3.5]. 
Figs. 13(b,f) show that near extrema, the MG3 constraint always provides “plenty of 
room” so that it has no effect. Note that for the three mesh intervals around the 
maximum, the interpolants in Figs. 13(b,f) are respectively identical to those of the 
unconstrained methods in Figs. 13(a,e). However, Figs. 13(d,h) show that the MS3 
constraint occasionally does not provide “enough room”, especially for a coarse mesh. 
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Figs. 10. Interpolation curves for FYitsck- Carlson’s data. 
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Figs. 12. Interpolation curves for f(x ) = e : 
on the domain [—2.8, 3.6], type I methods. 


52 







(fl) MP-quartlo 


o-t 


— 3 . 


7 . 


u . t i r 

-*• ttO 


MS3 — queurtlo 


4. 


Figs. 13. Interpolation curves for f(x) = e 
on the domain [-2.9, 3.5], type I methods. 
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(h) M3— B 
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Although the interpolants are “somewhat flat” at the local maximum in these figures, 
they do reach a strict local maximum at an interior point of the mesh interval. These 
results are consistent with the analysis in subsection 3.5. Also note that the MP 
constraint produces even “flatter” interpolants in Figs. 12(c,g), and the local maximum 
of the interpolants occur at a mesh point. 

In Figs. 14, the data are the same as those of Figs. 13; but the methods presented 
are of type II: the second and third-order limiter methods listed above in Figs. 9 and 
11. Since the data are smooth, the interpolants of the M3 methods in the right column 
are considerably smoother than those of the M2 limiters in the left. Note that the 
interpolants given by the “Superbee” and the M3-B methods in Figs. 14(g) and (h) are 
a little “fuller” than the others in the same column. As shown below, they also have 
larger errors, that is, the “visually pleasing” properties of the interpolants by these 
methods sire obtained at the cost of losing some accuracy. 

In the next examples, for comparison, the domain is identical to that of Hyman: 
[-1.7, 1.9]. In Table 1, the root mean square error, 


RMS error 




of the interpolants is compared as the mesh is refined from N = 8 to N = 64. The last 
column shows the result for an irregular mesh with N = 32, Aar m+3 /2 = 3Ax m+ i/ 2 , 
and every other mesh interval is of the same length. For type I methods, the MG 3 
constraint has no effect, while the MS3 constraint has no effect only when the mesh 
is fine enough, that is, when N > 16. It can be seen that the FD4 formula is highly 
accurate for a uniform mesh; however, for an irregular mesh, it is much worse than 
the parabolic and the quartic formulas. Similarly, the average formula is identical to 
the MP- parabolic method for a uniform mesh, but when the mesh is irregular, it has a 
larger error. Observe that the MS3, MG 3-parabolic, and all the M3 methods are more 
accurate than the M2-limiter methods. Among the limiters, those with g (1) 1/2 aie 

more accurate than those without. 

Note that for the minmod (or “Superbee”) limiter, /, — f'(xi ) and /,+ 1 — f'(xi+ 1) 
are of the same sign in the monotone part of the data since each fi is closer to (or further 
from) 0 relative to the exact solution f'(xi)- Consequently, Pf(x ) — f(x) changes sign 
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in [xj,x,-+i], and the root mean square error favors these methods. A more accurate 
norm is the average error of /,•: 

Average error = — — | /, - /'(*.•)!• 

t=m 


Method 

N = 8 

N = 16 

N = 32 

N = 64 

N = 32 irr. 

parabolic 

4.2E-3 

4.1E-4 

4.3E-5 

4.9E-6 

5.0E-5 

MG3-parabolic 

4.2E-3 

4.1E-4 

4.3E-5 

4.9E-6 

5.0E-5 

MS3-paxabolic 

4.9E-3 

4.1E-4 

4.3E-5 

4.9E-6 

5.0E-5 

MP-parabolic 

5.9E-3 

2.3E-3 

8.0E-5 

2.5E-5 

8.2E-4 

quartic 

3.4E-3 

7.4E-5 

2.3E-6 

1.2E-7 

7.5E-6 

FD4 

3.4E-3 

7.4E-5 

2.3E-6 

1.2E-7 

2.0E-3 

MG3-FD4 

3.4E-3 

7.4E-5 

2.3E-6 

1.2E-7 

2.0E-3 

MS3-FD4 

4.5E-3 

7.4E-5 

2.3E-6 

1.2E-7 

2.0E-3 

MP-FD4 

5.6E-3 

2.4E-3 

6.9E-5 

2.4E-5 

2.0E-3 

M2-Limiters: 

Minmod 

1.1E-2 

2.9E-3 

5.3E-4 

1.2E-4 

1.3E-3 

Van Albada 

9.9E-3 

2.6E-3 

3.3E-4 

6.5E-5 

1.2E-3 

Average 

5.9E-3 

2.3E-3 

8.0E-5 

2.5E-5 

1.1E-3 

Superbee 

7.9E-3 

3.3E-3 

4.7E-4 

1.2E-4 

1.2E-3 

Ftit sch-B ut land 

6.8E-3 

2.4E-3 

1.5E-4 

4.5E-5 

1.1E-3 

M3- 

MD or UNO 

5.2E-3 

4.4E-4 

4.8E-5 

5.3E-6 

6.2E-5 

V 

4.5E-3 

3.5E-4 

3.7E-5 

3.4E-6 

8.1E-5 

A 

4.3E-3 

4.1E-4 

3.8E-5 

3.4E-6 

8.5E-5 

B 

5.9E-3 

8.6E-4 

9.2E-5 

1.0E-5 

1.7E-4 

M4- 

MD 

3.2E-3 

1.7E-4 

8.1E-6 

4.4E-7 

1.3E-5 

V 

2.4E-3 

1.2E-4 

3.9E-6 

1.8E-7 

1.0E-5 

A 

2.7E-3 

1.8E-4 

5.5E-6 

1.8E-7 

1.4E-5 

B 

3.4E-3 

2.6E-4 

9.6E-6 

5.1E-7 

2.1E-5 


Table 1: RMS errors for f(x ) = e ** 
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Method 

N = 8 

N = 16 

N = 32 

N = 64 

N = 32 irr. 

parabolic 

5.8E-2 

1.6E-2 

4.0E-3 

9.9E-4 

3.1E-3 

MG3-paxabolic 

5.8E-2 

1.6E-2 

4.0E-3 

9.9E-4 

3.1E-3 

quartic 

4.1E-2 

1.7E-3 

9.9E-5 

6.1E-6 

7.4E-5 

M4-quartic 

4.8E-2 

1.7E-3 

9.9E-5 

6.1E-6 

7.4E-5 

M3-quartic 

5.4E-2 

1.7E-3 

9.9E-5 

6.1E-6 

7.4E-5 

MP-quartic 

6.2E-2 

2.3E-2 

8.6E-4 

3.9E-4 

7.4E-3 

M2-Limiters: _ . 

Minmod 

1.8E-1 

1.0E-1 

4.7E-2 

2.4E-2 

5.0E-2 

Van Albada 

1.3E-1 

6.0E-2 

1.7E-2 

5.7E-3 

2.8E-2 

Average 

7.6E-2 

3.6E-2 

4.7E-3 

1.4E-3 

2.8E-2 

Average Rational 

9.8E-2 

4.3E-2 

9.0E-3 

2.9E-3 

2.7E-2 

Average Cubic 

7.4E-2 

4.2E-2 

1.0E-2 

4.7E-3 

2.8E-2 

Harmonic Mean 

1.1E-1 

4.7E-2 

1.1E-2 

3.7E-3 

2.7E-2 

Fritsch-Butland 

7.6E-2 

3.5E-2 

9.0E-3 

5.8E-3 

2.8E-2 

Superbee 

9.9E-2 

7.7E-2 

4.0E-2 

2.2E-2 

4.7E-2 

M3- 

MD or UNO 

7.1E-2 

1.8E-2 

4.6E-3 

1.1E-3 

3.7E-3 

V 

5.5E-2 

9.5E-3 

2.5E-3 

5.6E-4 

2.8E-3 

A 

5.6E-2 

1.1E-2 

2.6E-3 

5.6E-4 

3.0E-3 

AR 

5.3E-2 

1.1E-2 

2.6E-3 

5.6E-4 

2.9E-3 

AQ 

6.3E-2 

1.4E-2 

2.9E-3 

5.8E-4 

3.2E-3 

HM 

5.2E-2 

1.0E-2 

2.5E-3 

5.6E-4 

2.9E-3 

FB 

6.0E-2 

1.5E-2 

3.9E-3 

9.0E-4 

4.2E-3 

B 

7.9E-2 

2.8E-2 

7.3E-3 

1.8E-3 

7.3E-3 

M4- 

MD 

4.2E-2 

5.3E-3 

6.1E-4 

7.6E-5 

4.6E-4 

V 

2.8E-2 

2.4E-3 

1.5E-4 

1.2E-5 

1.8E-4 

A 

3.4E-2 

3.5E-3 

1.9E-4 

1.2E-5 

2.4E-4 

AR 

3.3E-2 

3.5E-3 

1.9E-4 

1.2E-5 

2.4E-4 

AQ 

3.9E-2 

3.6E-3 

1.9E-4 

1.2E-5 

2.4E-4 

HM 

3.2E-2 

3.5E-3 

1.9E-4 

1.2E-5 

2.4E-4 

FB 

3.8E-2 

4.4E-3 

3.3E-4 

3.3E-5 

3.4E-4 

B 

4.8E-2 

6.9E-3 

7.2E-4 

8.5E-5 

6.5E-4 


Table 2: Average errors of fi for f(x) — e 
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In Table 2, the average errors are listed. As the mesh spacing is reduced by a factor 
of 2, the error of the parabolic method is reduced by a factor of roughly 4, the quartic 
method by 16, the minmod and “Superbee” limiter by 2, and the M3 methods by 4. 
The error of the M3-A method is about half that of the UNO, M3-B, and parabolic 
methods. This result is consistent with the analysis in subsection 3.5. Also note that 
limiters with p'(l) = 1/2 are more accurate than those without. Among limiters with 
this property, the average errors are essentially the same. The errors of the M4-MD 
and M4-B methods reduce by a factor of roughly 8 as the mesh is refined, and for the 
M4-A method, the factor is larger than 8 in this case. The M3-B and M4-B methods 
have the largest error in the M3 and M4 classes, respectively. 

The examples in this subsection show that accurate, monotonicity-preserving, and 
“visually pleasing” interpolants are produced for smooth data by the following methods: 
MS3, MG3, M3-parabolic and quartic, M3-V, M3-A, M4-V, M4-A, and M4-quartic. 
These methods, except for the M3-V and M4-V, also produce monotonicity-preserving 
and “visually pleasing” interpolants for rough data, as shown in the previous subsection. 

7. Discussion and conclusions. The above methods can be applied to con- 
servation laws, see [16], [29] for some preliminary results in the piecewise linear case. 
We hope to report further results for the piecewise linear and parabolic methods in a 
forthcoming paper. 

For rough, smooth, monotone or nonmonotone data sets, the methods discussed 
produce stable, monotonicity-preserving, and “visually pleasing” C 1 Hermite cubic 
interpolants which are at least third-order accurate. The two simplest methods are 
the MS3-parabolic and M3-A methods, where the latter has an average error for the 
derivatives roughly half that of the former. Without increasing the stencil, the M3- 
quartic method is fourth-order accurate— the highest possible order of accuracy for 
a cubic interpolant — except adjacent to strict local extrema and inflection points, 
where the accuracy degenerates to third-order. Uniform fourth-order accuracy can be 
obtained by using the M4-A or M4-quartic methods. 
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